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ABSTRACT 

A massive black hole (MBH) consumes stars whose orbits evolve into the small phase-space volume of un¬ 
stable orbits, the “loss-cone”, which take them directly into the MBH, or close enough to interact strongly with 
it. The resulting phenomena: tidal heating and tidal disruption, binary capture and hyper-velocity star ejection, 
gravitational wave (GW) emission by inspiraling compact remnants, or hydrodynamical interactions with an 
accretion disk, are of interest as they can produce observable signatures and thereby reveal the existence of the 
MBH, affect its mass and spin evolution, probe strong gravity, and provide information on stars and gas near 
the MBH. The continuous loss of stars and the processes that resupply them shape the central stellar distri¬ 
bution. We investigate relativistic stellar dynamics near the loss-cone of a non-spinning MBH in steady-state 
analytically and by Monte Carlo simulations of the diffusion of the orbital parameters. These take into account 
Newtonian mass precession due to enclosed stellar mass, in-plane precession due to general relativity, dissi¬ 
pation by GW, uncorrelated two-body relaxation, correlated resonant relaxation (RR) and adiabatic invariance 
due to secular precession, using a rigorously derived description of correlated post-Newtonian dynamics in the 
diffusion limit. We argue that general maximal entropy considerations strongly constrain orbital diffusion in 
steady-state, irrespective of the relaxation mechanism. We identify the exact phase-space separatrix between 
plunges and inspirals, predict their steady-state rates, and verify they are robust under a wide range of assump¬ 
tions. We derive the dependence of the rates on the mass of the MBH, show that the contribution of RR is 
small, and discuss special cases where unquenched RR in restricted volumes of phase-space may affect the 
steady-state substantially. 

Subject headings: black hole physics — galaxies: nuclei — stellar dynamics 


1. INTRODUCTION 
1.1. Background 

Strong interactions of stars with a massive black hole 
(MBH) in a galactic center lead to a variety of extreme phe¬ 
nomena, as well as provide mass for the growth and evolution 
of the MBH. The small phase-space volume of orbits whose 
periapse lies close enough to the MBH to lead to a strong in¬ 
teraction is called the loss-cone prank & Rees||l976| , since 
in most cases the interaction destroys the star, either immedi¬ 
ately (for example by a direct plunge through the event hori¬ 
zon or by tidal disruption outside it (Rees] 1988j i or gradually, 
after the orbit decays by some dissipation mecha nism (for ex¬ 
amp le the emission of gravitational waves (GW) (Hils_& Ben- 
der|1995|l, tid al heating of the star by the MBH ( Alexander & 
Morris 3003), or drag against a massive accretion disk (|Os- 
triker|1 983)). Even when the star s urvives the encounter, fo r 
example in a tidal scattering event (Alexander & Livi oj200T| >, 
the restricted set of orbits that allow such near-misses lie very 
close to the loss-cone phase-space. Since typically stars do 
not survive long on loss-cone orbits, the key questions is how, 
and at what rate are these orbits repopulated by new stars. The 
stellar dynamical study of this question is known as loss-cone 
theory. A main re-population channel is by dynamical relax¬ 
ation mechanisms, which randomize stable orbits and causes 
them to diffuse in phase-space into the loss-con^] The close 
interaction event rates in the steady-state of dynamically re¬ 
laxed systems are of particular interest, both because these 
can be derived from first principles independently of initial 
conditions, and because these correspond, statistically, to the 
cases most likely to be observed. 


1 Other possibilities include e.g., in-situ star formation or galaxy mergers. 


Studies of the loss-cone can be broadly categorized by four 
criteria: whether they deal with processes that lead to im¬ 
mediate stellar destruction (infall) or a gradual one (inspi¬ 
ral); whether they are strictly Newtonian or include also gen¬ 
eral relativity (GR), fully or perturbatively; whether they in¬ 
clude only slow non-coherent two-body relaxation (i.e., non¬ 
resonant relaxation, NR |Chandrasekhar|1944[ > or also fast co- 
herent relaxatio n (known as resonant relaxation, RR, R auch &| 
Tremaine|1996| See Section^; and finally by the calculation 
methods employed: whetherthey are analytical, or based on 


the diffusion approximation either by direct numerical solu¬ 
tions of the Fokker-Planck (FP) equations or by Monte-Carlo 
(MC) methods, or whether they employ direct iV-body simu¬ 
lations. 

Early studies focused on the infall rates of tidal disruption 
events in the Newtonian approximation, using analytic and 


FP-based methods ((Frank & Rees|| 1976; 

Young et al. 1977; 

Cohn & Kulsrud 1978, Shapiro & Marc 

hant 1978). These 

studies were subsequent] 

y updated and generalized to include 


some deviations fro m sp herical symmetry (Syer & Ulmer 
Magorrian & Tremaine|1999[ see review by Alexander 


1999 


2012). The prospects of detecting low-frequency GW from 


compact remnants spiraling into MBHs (extreme mass ratio 
inspiral (EMRI) events) with long-baseline space-borne GW 
detectors ( |Amaro-Seoane et al.|2007 1 motivated studies of the 
inspiral rates for EMRIs in using FP-based methods in the 


NR-only limit with perturbative GR (Hils & Bender 

1995 

Sigurdsson & Rees 1997; Miralda-Escude & Gould 

2000 

Freitag 2001 2003 [[Ivanov 2002; Hopman & Alexander 

2005 

2006b[ see review by Sigurdsson 2003 ). The ettects ot the 


Seoane et al. |2013| l. 

A unifying framework relating plunge and inspiral pro- 
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cesses was formulated by Alexander & Hopman) (J2003j) and 

G'c 


used to estimate infall and inspiral event rates in the Galac- 
tic Center in the NR-only limit: infall by direct plunge and 
tidal disruption, inspiral by GW emission and tidal heat¬ 
ing ( Alexander & Morris||2003|), as well as tidal scattering 
events ( {Alexander & Livioj 25011. Different attempts to es¬ 
timate the infall and inspiral rates yielded a wide, uncertain 


range of values that spans several orders of magnitudes (Sig- 
|urdsson|2003| |Aiexander|2012| l. 

Fast relaxation by RR can be effective on the small spatial 
scales where most EMRIs originate and importantly, where 
stellar orbits are observed in the Galactic Center and can 


therefore provide empirical constraints (Hopman & Alexan- 
|der|2006a| . This realization motivated a re-evaluation of re¬ 
laxation processes and their impact on dynamics very close 
to MBHs. An approximate comparison of the relative rates 
of RR and NR suggested that the branching ratio between 
plunges and inspirals depends strongly on the efficiency of 
RR, which was then poorly u nderstood (|Hopman & Alexan- 


der|200 6a Eilo n et al.|2009 1. This added a yet larger uncer¬ 


tainty to EMRI rate estimates. A key question is the physical 
origin and characteristics of the quenching mechanism that 
perturbs the near-Keplerian symmetry that generates RR, and 
causes the orbits to drift in phase-space from their initial val¬ 
ues. 


Initial analysis of RR in the relativistic context (Rauch & 
Tremaine][T996 ) indicated that rapid GR precession on very 


eccentric orbits likely plays a key role in quenching RR. Im¬ 
portantly, GR quenching can prevent RR from rapidly pushing 
all stars into plunge orbits, thereby allowing slow inspiral to 
produce detectable p eriodic GW signals (EMRIs) (Hopman 


|& Alexander|2006a) l. In these earlier studies the determinrstrc 
GR precession was treated as an effective stochastic perturba¬ 
tion of the Keplerian orbits. 

First indications that the precession cannot be treated that 
way, and that the long-timescale behavior of RR is not well 
described as a Markov process (random walk), were un¬ 
covered in post-Newtonian small IV-body simulations of di¬ 


rect plunge and GW inspiral events (Merritt et al. 2011 
MAMW11). These revealed oscillatory orbital behavior at 
high eccentricities that appeared to act as a barrier against 
further evolution to even higher eccentricities and subsequent 
infall or inspiral. MAMW11 dubbed this dynamical phe¬ 
nomenon the Schwarzschild Barrier (SB), and showed that the 
oscillations are well approximated by the simple ansatz of as¬ 
suming 2-body GR dynamics in the presence of a randomly 
oriented fixed force vector ( |Alexander||2010) >, representing 
the residual force due to the background stars. While the 
effect appeared related to the EMRI-preserving RR quench¬ 
ing predicted by Hopman & Alexander (2006a). its magni¬ 
tude seemed much stronger than anticipated, in that it not 
only damped the RR torques, but actually appeared to pre¬ 
vent the orbits from interacting closely with the MBH_at all. 


The larger-scale relativistic IV-body simulations of Brem et al. 
(20141 BAS 14) confirmed that GR precession quenches RR 


roughly on the scale of the SB, and concluded that the result¬ 
ing EMRI rates are consistent with those driven by NR. 

The SB phenomenon was subsequently explained rigor¬ 
ously in terms of the adiabatic invariance (Al) of the GR pre¬ 
cession against the coherent RR torques, when the precession 
period is shorter than the typical RR coherence time (the rj- 
formalism,(Bar-Or & AlexanderJ2014i see review by Alexan¬ 
der 2015|. By describing the RR torques due to the back¬ 


ground stars in terms of a correlated noise field, it is possi 


ble to formulate an effective FP description for RR that takes 
into account Al, and to derive the corresponding effective 
diffusion coefficients (DCs), whose form and behavior de¬ 
pends critically on the assumed temporal smoothness of the 
noise model. The continuous orbital evolution of the stellar 
background suggests that the physically correct form of the 
stochastic torques is that of a smooth (infinitely differentiable, 
C °°) noise. The Al is maximal for smooth noise. In that case 
its dynamical effect can be described as a faster than expo¬ 
nential suppression of the diffusion coefficients below some 
critical angular momentum limit. The vanishing phase-space 
density past this limit grows so slowly (~ log i), that the limit 
can be considered as an effective barrier fixed in time. While 
this limit is not a true barrier, nor a reflecting one, it does 
effectively divide phase-space into a region where RR can be 
efficient, and a region where it cannot. As we show below, the 
presence of the competing process of NR substantially limits 
the significance of Al in the dynamics of the loss-cone on long 
timescales (of order of the NR relaxation time). 

This study focuses mainly on the implications of NR and 
RR around a MBH for loss-rates. However, these dynamical 
processes are also relevant for understanding and modeling 
other processes around MBHs, and in particular the Galactic 
MBH, SgrA*. Although the inner Galactic Center contains a 
relatively small and manageable number of stars by the stan¬ 
dards of current Newtonian IV-body codes, it is still extremely 
challenging to simulate it directly both because of the extreme 
dynamical range introduced by the high MBH to star mass ra¬ 
tio, and because of the added complexity of the GR equations 
of motion. The impact of MBH spin and RR on orbital tests of 
GR in the Galactic Center were studied with post-Newtonian 
small iV-body simulations ( Merrit fet al.||20IQ| . A study of 
the implications of RR for the formation mechanisms of the 
of stars orbiting SgrA*, either resorted to large Newtonian- 
only IV-body simulations (Perets et al . [2008) 1, or substantially 
underestimated the efficiency of Al in quenching RR by us¬ 
ing a MC scheme based on the simple fixed force ansatz with 
non-differentiable (C°) noise, to study the implications of the 
SB for stars in the Galactic Center ( jAntonini & Merritt|2013[ 
|Antonini|2014| l. 


1.2. Objectives and oven’iew 

The objectives of this study are to integrate the recent in¬ 
sights about the role correlated noise plays in determining the 
properties of RR and its formulation as an effective diffusion 
process ( |Bar-Or & Alexander|2014| > together with the known 
properties of NR; to derive a rigorous computational frame¬ 
work for calculating the steady-state phase-space density near 
the relativistic loss-cone and the resulting loss-rates, and to 
use this framework for a systematic study of the dependence 
of the results on the various physical mechanisms involved in 
the dynamics: mass (Newtonian) precession, GR precession, 
GW dissipation, the RR noise model and coherence time. The 
ultimate objective is to provide well-defined estimates of the 
infall and inspiral rates (including, but not limited to direct 
plunges, tidal disruptions and EMRIs) and their scaling with 
the properties of the galactic nucleus (MBH mass and stel¬ 
lar density). These can then inform design decisions about 
planned surveys and experiments, and serve as benchmarks 
for more detailed future studies. 

We focus our study on a simplified galactic nucleus con¬ 
taining a stationary non-spinning (Schwarzschild) MBH sur¬ 
rounded by a Keplerian, spherically symmetric (in the time- 
averaged sense), power-law cusp of single mass stars (the 
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background cusp). Direct relativistic iV-body simulations 
generate, by construction, the correct dynamics, but are 
presently limited by computational costs to unrealistically 
small N, which generally cannot be scaled up to astrophys- 
ically relevant values since different dynamical mechanisms 
scale differently with N (e.g., Heggie & Hut [20031 . More¬ 
over, they allow little freedom to switch on/off the various 
physical mechanisms that affect the outcome. It is therefore 
difficult to disentangle their contributions and interpret the re¬ 
sults. 

Here we follow a different approach. We represent the evo¬ 
lution of the system in the realistic large-iV limit as a super¬ 
position of diffusion processes. This enables us to isolate and 
study the effect of the different dynamical mechanisms, and 
thereby obtain an analytic description of the system. 

We calculate the loss-cone phase-space density and loss- 
rates by two complementary methods. We show that the dif¬ 
fusion in phase-space is well approximated as a separable 
process: fast diffusion in angular momentum, superposed on 
a slow diffusion in energy. We then use this separation of 
timescales to derive analytically the steady-state properties of 
the system. We also solve the diffusion in energy and angu¬ 
lar momentum phase-space numerically by MC simulations, 
which are statistics-limited, but have the advantage of flex¬ 
ibility in introducing additional dynamical effects and con¬ 
straints. We cross-validate these two calculation methods, 
and also compare the MC results to the A'-body loss-rates of 
MAMW11 and BAS 14, and reproduce the AI effects of |Bar-| 
|Or & Alexander| ( f2014 1 in the absence of NR. 

This paper is organized as follows. In Section [2] we re¬ 
view the dynamical process of relaxation near a MBH in a 
galactic nucleus. We present a unified framework for describ¬ 
ing both non-coherent two-body relaxation and coherent res¬ 
onant relaxation. We discuss the role of secular processes in 
the emergence of adiabatic invariance in the long-term orbital 
evolution of the system. In Section[3]we describe the structure 
and properties of phase-space near the loss-cone, and derive 
analytic estimates for the stead y-state distribution and loss- 
rates. We start, in Section by formulating the diffusion 
equ ation s which govern the evolution of the system. In Sec¬ 
tion |3.2| we describe of the diffusion process in terms of the 
streamlines of the probability flow, which provide a power¬ 
ful visua l represe ntation of the dynamics, and guides us in 
Sections |3.3f|3.4| in solving the steady-state distribution and 
loss-rat es (a t this stage, without RR or GW dissipation). In 
Section [X5j we show that GW dissipation separates the prob¬ 
ability flow into two distinct regions in phase-space: a region 
where stars can inspiral into the MBH while emitting GWs, 
and a region where stars plunge directly into it. The inspiral 
event rate is then calculated exactly by locating the separa- 
trix that d emar cates the two regions (Appendix [A]). Finally, 
in Section [376] we show that RR has only a smalTimpact on 
the steady-state density and loss-rates, and provide a method 
to quantify its effect. In Section [4] we present our MC proce¬ 
dure for modeling orbital evolution in phase-space (The MC 
procedure and the NR DCs th at en ter it are described in Ap¬ 
pendices IB] and |Cj. In Section 4.1 we validate the description 
RR and the emergence of AI in angular momentum-only sim¬ 
ulations against the analytic results of Bar-Or & Alexander 
(j2014j), and show that AI is very efficiently suppressed by NR 
on long timescales. We also show that over short timescales, 
AI induces the “Schwarzschild Barrier” phenomenon seen in 
angul ar momentum and energy phase-space (Merritt et al. 


|2011)>, and demonstrate that this dynamical feature is erased 


over long timescales. We compare in Section|4.2|the MC code 
against the small A'-body loss-rates of Merritt et al.| (|2011 1 
and |Brem et al.| (|2014)», and show that the MC and Tv -body 
give consistent results (The derivation of steady-state rate esti¬ 
mates from MC simulations is described in Appendix |B.3[ ). In 
Section [5] we explore the robustness of the MC-derived rates 
under v ariou s dynamical approximations and assumptions. In 
Section |5.1| we estimate the rates for the prototypica l tar get 


of low-frequency GW, the Galactic Center. In Section 5.2 


we 


derive analytically, and confirm with the MC simulations, the 
weak scaling of the rates with the MBH mass. We di scuss and 
summarize our results in Section [6] In Section 16.11 we focus 
on the role of the principle of maximum entropy (Appendix[E]i 
as a guiding principle in the derivation of the DCs. We argue 
that RR typically does not play a major role in the steady- 
state dynamics of the loss-cone. We illustrate this analysis by 
presenting a fine-tuned idealized counter-example where RR 
may substantially affect the loss-rates: the interaction of icy 
planetesimals with a massive circumnuclear accretion disk. 
We con clud e by discussing the limitations of our analysis in 
Section [C2l 


2. RELAXATION AROUND A MASSIVE BLACK HOLE 

We present here an overview of dynamical relaxation 
around a central MBH, and in particular of resonant re¬ 
laxation, using two complementary approaches. The first, 
which is closer to the conventional description of the sub¬ 
ject (e .g., |Rauch & Tremaine 1 1 996[ [Hopman & Alexa nder! 
2006a}, serves to introduce the basic terms and ideas, and con¬ 
nect the present work to past studies. The second presents a 
unified framework for describing the coherent and stochastic 
aspects of relaxation, using adiabatic invariance as the uni¬ 
fying concept in all forms of relaxation. This connects the 
present work to recent results on the representation of RR as 
an effective diffusion process, a key tool used here to investi¬ 
gate loss-cone dynamics in the large-A limit. 


2.1. Two-body relaxation 

To understand the relation between NR and RR, consider 
a spherical stellar system composed of stars of mass M* or¬ 
biting a central massive object M . A> M„, and focus on a 
test star at radius r from the center, where the local stellar 
number density is n* (r). The average net force exerted on 
the test star by the diV* (< b ) ~ n*b 2 db background stars 
in a thin small shell around it with radius b -C r and width 
d b <C b, is zero. However, the Poisson fluctuations in the 
positions of these d./V* discrete masses generate a residual 
force of magnitude yj ( F 2 ) ~ ^/dN+GM 2 /b 2 . This force 
persists in direction and magnitude until the stars generating 
it move substantially. For a random stellar velocity field with 
dispersion er 2 ~ GM t /r, this coherence time is T^ R ~ 
b/o <A r/a ~ P, where P is the orbital period. Since 
T^ r <C P, the net encounter is impulsive—a collision, and 
so in the case of NR, the coherence time is the collision time. 
The change in velocity due to the residual force over time 
T^ r is Sv ~ yj(F 2 )T^ r /M*.. Over times t > T^ R , these 
impulses add non-coherently, and the accumulated change in 
velocity per unit time is (An 2 ) t ~ (G 2 M 2 n ic /a) db/b. In¬ 
tegration over all shells from to 6 rnax < r (assuming 
n* is constant) yields the NR diffusion timescale 7’v u ~ 
o- 2 /(Av 2 ) t ~ a 3 / (G 2 M 2 n* log A) where A = 6 max /6 min 
is the Coulomb factor. These local changes in the test star’s 
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velocity lead to changes in orbital energy and angular momen¬ 
tum at a rate (A E 2 ) /E 2 ~ (AJ 2 ) / J 2 ~ I/T^r, where 
J c = VGM^ is the circular angular momentum, and a the 
semi-major axis (sma). 


2.2. Resonant Relaxation 

In a nearly-symmetric potential where the background or¬ 
bits are nearly fixed on timescales T c P , the test star 

interacts over a long period of time with the entire phase- 
averaged background orbits, and not only instantaneously 
with the segment of the orbit closest to it. That is, the in¬ 
teraction is non-local. As in the NR case, the discrete num¬ 
ber of stars on the scale of the test star’s orbit, TV* (< r), 
gives rise to random fluctuations in the force on it, \J(E 2 ) ~ 
yj TV* (< r)GM 2 /r 2 , which persist on timescale T c . How¬ 
ever, unlike in the NR case, here T c P, and so the orbital 
energy is adiabatically conserved (see below) and the residual 
force affects only the orbital angular momentum. 


2.3. A unified description of relaxation 

Conventionally, orbital evolution is described by a com¬ 
bination of two distinct classes of processes: stochastic re¬ 
laxation processes (two-body relaxation) and coherent (sec¬ 
ular) processes (RR). Here we present a unified framework 
(summarized in Table IT) for describing and analyzing all the 
dynamical processes that drive orbital evolution in terms of 
short-timescale coherent processes that effectively contribute 
as stochastic processes on longer timescales. 

We begin by considering the case where the test star is 
statistically indistinguishable from the background stars (this 
case was the focus of early works on RR, and is recast here in 
a more general context). 

It is instructive to consider the various relaxation processes 
in the order of their associated coherence times. When two 
stars interact impulsively (as in the case where the impact pa¬ 
rameter b is much smaller than the sma a of the orbit around 
the MBH), the interaction time is effectively limited to the 
crossing time of the closest approach. Since during the inter¬ 
action, the force on the test star is nearly constant, the duration 
of the interaction (the collision timescale) is the coherence 
timescale T^ R ~ b/o of two-body relaxation. As long as 
only interactions with small enough b < a are considered, so 
that T f MR < P(a), the Hamiltonian cannot be orbit-averaged, 
and therefore all the orbital elements can change by the inter¬ 
action. 

The treatment of two-body relaxation is based on the ap¬ 
proximation Tfi R —► 0, that is, that the interaction time is 
shorter than any other relevant timescale in the system, and 
that individual collisions are uncorrelated. In that limit the 
process is Markovian (Nelson & 'I remaine {1999)[B~ar-Or et ah] 
2013)1 and can therefore be described as diffusion in phase- 
space. 

Two-body interactions with a large impact parameter b > a 
(i.e., soft encounters) such that T^ R ~ b/a > P(a ) are no 
longer impulsive. This means they cannot be described as oc¬ 
curring instantaneously and locally between two point parti¬ 
cles, and therefore can no longer be described by the standard 
two-body relaxation formalism. In particular, the interaction 
is no longer Markovian, since the test star is repeatedly af¬ 
fected by the same perturbin g sta r. Since it can be shown (for 
energy relaxation, Bar-Or et al.|2013) that these soft two-body 
encounters do not contribute much to the total relaxation, an 


approximate cutoff on the maximal impact parameter is intro¬ 
duced via the Coulomb logarithm term. 

Two-body interactions in the extreme soft limit can be de- 
scribed in terms of an effective diffusion process ( |Bar-Or &| 
Alexander 2014). Since T c > P, the Hamiltonian can be 
double-averaged over both the orbit of the test star and the 
orbits of the background stars. The averaged Hamiltonian is 
then independent of the mean anomaly, and so the Keplerian 
energy (semi-major axis) is adiabatically conserved—proof 
that the contribution of soft collisions to energy relaxation is 
negligible. The double-averaged Hamiltonian no longer de¬ 
scribes point particles, but rather interaction between Keple¬ 
rian ellipses (“mass wires”), which mutually torque each other 
and exchange angular momentum, but not energy. Here, the 
force on a test ellipse by the background ellipses remains con¬ 
stant as long as the orbital orientations of the background el¬ 
lipses remain fixed (i.e., over the coherence time T RR ). In 
analogy to the case of point-point two-body relaxation, this 
coherence time can be considered as the interaction (“colli¬ 
sion”) time. The coherence time is determined by the fastest 
process that can reshuffle the background orbital orientations. 
Since typically the background stars are not on relativistic or¬ 
bits, the dominant shuffling process is the retrograde in-plane 
drift of the argument of periapse, ui, due the enclosed stel¬ 
lar mass inside the orbits, on the mass-p recession timescale 
Tm ~ QP/N (e.g., Hopman & Alexander 2006a i. As long 
as there are no competing processes with timescale shorter 
than Tm that could randomize the residual forces of the orbit- 
orbit interactions, these will torque the orbits and change their 
angular momentum in a coherent (oc t) fashion. Therefore, 
the Markovian assumption can used and the diffusion rate is 
(t 2 ) RR Tm ~ QP- This regime of RR is sometimes called 
“scalar RR” since it can change the magnitude of the angular 
momentum as well as its direction. 

On timescale substantially longer than the precession pe¬ 
riod, the Hamiltonian can also be double average over the 
argument of periapse and the interaction is then between 
mass annuli (K ocsis & Tremaine] |2015| ). In this case the 
collision (coherence) time is the self-decoherencing (or self¬ 
quenching) time Tf RR = T sq ~ J 2 / ( t2 ) ~ QP/s/N on 
which the annuli are re-shuffled by their own mutual torques. 
The residual torque (t 2 ) rr is now also averaged over the 
argument of periapse, which leads some cancellation of the 
torques, and therefore ( t2 ) vRR < ( T ‘ 2 ) RR an d the diffusion 
rate is (t 2 ) vRR Tf RR ~ QP/y/N. This regime of RR is 
sometimes called “vector RR” since it can change only the di¬ 
rection of the angular momentum, but not its magnitude. In 
general, as we consider longer coherence timescales, it be¬ 
comes possible to average over yet more degrees of freedom 
(phases), and the averaging results in an effective potential 
that is yet more symmetric. This in turn reduces the magni¬ 
tude of the residual forces. It is found that the countervailing 
effects of smaller torques but longer coherence times lead to 
faster relaxation T r oc v 2 / F 2 ) T c ), see TablejTj 

We now turn to the case where the timescales of the test- 
star and the background are different, which was not treated 
rigorously until recently ( Bar-Or & Alexander|2014 1 , follow - 
ing the discovery in V-body simulations (Merritt et al. |2011) 
of an abrupt transition in phase-space to adifferent dynami¬ 
cal regime where orbital evolution is governed by determin¬ 
istic rather than stochastic processes. In this case the pre¬ 
cession of the test star is due to the combined prograde GR 
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in-plane precession with period T G r and the retrograde New¬ 
tonian precession due to the stellar mass enclosed inside the 
orbit. Since the GR precession rate Tq^ diverges as 1 /j 2 , 
where j = \/"P e 2 and e is the eccentricity, eccentric stars 
with j smaller than some critical value jo will precess much 
faster than the background, i.e., T G r < T c . In this case, the 
Hamiltonian can be averaged over both the mean anomaly of 
the background and of test star, as well as over u> of the test- 
star. As a result, the test star’s j is adiabatically conserved. 
The transition between unconstrained RR-driven diffusion at 
j jo to strict adiabatic invariance at j < jo was calculated 
for several models of the effective background “noise” (resid¬ 
ual torques) in terms of effective diffusion coefficients by Bar 


Or & Alexander (2014 1 , who showed that for a smoothly vary¬ 
ing background the transition at jo is extremely sharp, with 
the j-diffusion coefficients suppressed exponentially with the 
argument (T c /T G r) 2 cx (jo/j) 4 . 


3. THE PHASE SPACE OF THE LOSS-CONE 

We describe the dynamics of the loss-cone in the Keplerian 
approximation, where the gravitational potential of the stars is 
assumed to be negligible relative to that of the central MBH. 
In that limit, the Keplerian orbital energy is E = GM, /2a, 
where the stellar dynamical convention E = —E trU e > 0 
for bound orbits is adopted. The orbital angular momentum 
is parameterized by j = J/J c (a). We marginalize the dy¬ 
namics over the orbital angles and consider the evolution in 
the (a, j) phase-space. We assumes a stationary non-spinning 
MBH of mass M, that is surrounded by an isotropic power- 
law cusp of stars of mass M* each with number density profile 
n* = no{r/ro)~ a . The mass ratio is denoted Q = M, / M„. 
The cusp is assumed to extend between a,,,;,, and a max . The 
inner boundary at a m j n is an absorbing boundary, set at the 
innermost stable circular orbit. The outer boundary at a max is 
the interface to the galaxy around the central cusp, which pro¬ 
vides an effectively infinite reservoir of stars to replace those 
that are lost into the MBH, or evaporate back to the galaxy. 
The Newtonian gravitational dynamics are described in terms 
of DCs, whose functional form and normalization are calcu¬ 
lated assuming the background cusp and an isotropic distribu¬ 
tion of angular momentum. The DCs, together with additional 
non-Newtonian processes such as an absorbing boundary at 
the last stable orbit (LSO) and GW dissipation of energy and 
angular momentum, are then used to generate the dynamics 
of test particles, and derive their steady-state loss-rates and 
phase-space density. 

Stars reach the MBH by crossing the LSO loss-line in 
(a, j) phase-space (Figure lk at ji c (a ) = 4 /\Ja/r g , where 
r g = GM., /c 2 (This value of Ji c = J c ji c is exact for a zero 
energy orbit). In addition, it is useful to define in the statistical 
sense the locus of “no-return” for GW inspiral (EMRI). Con¬ 
ventionally, this is defined by a comparison of timescales as 
the locus where the time to spiral into the MBH by the emis¬ 
sion of GWs, t.Gw(a,j), is shorter than the time needed to 
scatter across the LSO line by NR, (j — ji c ) 2 Tj(a, j), where 
Tj is the J-diffusion timescale (Appendix (A). Note that the 
GW timescale line shown in Figure [I] is calculated both us¬ 
ing the common simplification taw = j 2 Tj, and with the 
more accurate form t G w = (j — jic) 2 Tj , which yields an 
arc-like shape that peaks well below the point where the ap¬ 
proximate power-law GW line intersects the LSO line. The 
maximal sma along these lines, a G w, is then interpreted as 
the critical sma for EMRIs, below which phase-space trajec- 



FlG. 1.— The (a, j) phase-space of the loss-with the various critical lines 
and regions, for a model of the Milky Way model with Q = 4 X 10L mass- 
precession coherence time and a Gaussian noise model (Section |5.If Or¬ 
bits in the gray area below the LSO line (red) are unstable and promptly 
plunge into the MBH event horizon (plunge track example in light red line). 
Where RR diffusion is faster than NR diffusion (yellow region), RR domi¬ 
nates the dynam ics. The S-stars observed near the MBH of the MW (red cir¬ 
cles) (Gillessen et al.|2009b|i lie in the RR dominated region. AI suppresses 
RR torquing below the Al line (gray). Inside the phase-space region delim¬ 
ited by the GW line (blue), GW dissipation is faster than NR J-scattering 
and orbits spiral into the MBH by the emission of GW (inspiral track exam¬ 
ple in light blue line). The critical sma for EMRIs, a gw (thin black line), 
corresponds to the maximum of the GW curve; below it stars become EMRIs 
before they cross the LSO. The approximate power-law GW line with the 
often assumed simplification ji c —> 0 (dotted blue line), substantially over¬ 
estimates ciQ\y. The exact separatrix streamline (magenta) provides a more 
accurate estimate of a a w than either of the timescale-based GW lines. 

tories cannot (statistically) avoid crossing the GW line and be¬ 
coming E MRI s. The EMRI rate scales as cx a G w (Eq. ©). 
In Section [53] and AppendixlAlwe formulate a more rigorous 
criterion for the GW line by identifying the exact separatrix 
between phase-space streamlines that plunge directly into the 
MBH, and those that inspiral into it (see Figure [4}. This re¬ 
sults in an intermediate value of a G w- These three estimates 
of the GW line are plotted in Figure |T| for reference, and cor¬ 
respond to different EMRI rate predictions. It should be em¬ 
phasized that the GW line does not enter the MC procedure 
directly, but is an emergent property. Here we use the sep¬ 
aratrix method for analytic rate estimates, which accurately 
reproduce the MC results (Figure [T7). 

We derive here analytic estimatesfor the steady-state distri¬ 
bution and the flux of stars through the loss-cone, quantify the 
contribution of RR to the loss-rates, and validate our estimates 
by MC simulations. 


3.1. Diffusion equations 

On long-enough timescales, where relaxation can be de¬ 
scribed as a diffusion process (Bar-Or et al. 2013[ Bar-Or 
|& Ale xander 2014), the evolution of the probability density 
function, n (E, J, t ), in (E, J ) phase-space can be describe 
by an FP equation 


dn ( E , J, t) dS e ( E , J, t ) dSj (E , J, t ) 


dt 


dE 


dJ 


(1) 
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TABLE 1 

A UNIFIED FRAMEWORK FOR RELAXATION PROCESSES (SEE TEXT) 


Process 

Effective 

particles 

Averaged quantity 

Conserved 

quantities 

Coherence time 

T c 

Residual force magnitude 

Relaxation time 

V 2 / (T C F 2 ) 

NR 

Points 

None 

None 

t nr k p 

~ \fNLGMi, sJQ/ a 2 f\ 

~ Q 2 P/(N log Q) 

RR 

Ellipses 

Mean anomaly 

E 

P < Tfl R ~ T m < T p 

~ sjNLGM*/a 2 

~ QP 

vRR 

Annuli 

Argument of periapse 

E, J 

T p < Ty RR ~ T sq 

~ yjNLGMt/a 2 

~ qp/Vn 


“Integrated over all impact parameters 


where the probability current densities in the E and J “direc¬ 
tions” are 

Se {E, J , t) = D E n ( E , J,t) — - — [D EE n (E, J, t)] 

^2~dJ \P E jn(E,J,t)\ , (2) 

and 

Sj (E, J, f) = Djn (E, [ D JJ n ( E i J, *)] 

~2 We ^ Ejn ^ ,i )] ’ ( 3 ) 

and where Bb, .Dj-j, Dj and D E j are the DCs that 

describe the combined effect of NR and RR. 

Integrating over J from to J c we obtain the total proba¬ 
bility current density gradient (loss-rate) per unit energy 


dn (E, t) 

at 


/ J c f) Qj-, 

-gfdJ-Sj(E,J c )+Sj(E,Ji c ) , 

(4) 


in the presence of a loss term Sj (E, Ji c ), resulting from the 
flux of stars through the loss-cone, per unit energqq with the 
J-averaged diffusion coefficients 

D E (E)=n~ 1 (E) [ D e (E, J) n (E, J) dJ , (10) 

Jjlc 


D ee (E) = n~ 1 (E) [ ' D EE (E, J) n (E, J) dj . (11) 

JJlo 

3.2. Probability flow in phase space 

The effects of the various physical mechanisms are more 
clearly apparent in the flow patterns in phase-space. Since 
the physical flow is stochastic, it is more useful to describe 
it in terms of the flow of the probability density. In steady- 
state, the FP equation (Eq. {I}) can be written as a continuity 
equation of a compressible flow 

O r\ 

— [n (E, J) v E ] + -gj [n (E, J) vj\ = 0 , (12) 


where n {E) = J'j° n ( E , J) d.I is energy probability density. 

It is convenient to transform these expressions to (E,j) 
since for j = J / J c , the current density in the ^-direction is 
zero at the boundary (j = 1). Generally, under the coordi¬ 
nate change x —► x' (here (E, J) -A ( E,j )), the probability 
density currents, S, —> S[, transform as 




Bx 

Bx' 


M 

dx k 


Sk (x) . 


(5) 


Thus 

B J 

Sj (E,j) = Sj (E, J) -Jq^Se (E, J ) , (6) 

and since Sj ( E,j = 1) = 0, we have 


Bn (E, t ) 
dt 


-1 ° d ~§§dJ - H S E (E, J c ) + Sj (E, J lc ) , 

(7) 


and 


Bn (E, t) 
dt 


_B_ 

dE 


S E dJ + Sj (E, Jlc) 


( 8 ) 


Using Eq. d2| and the fact that D E j(E,J c ) = 

D ee (E, J c ) ffJjdE (Appendix Icb, allow us to obtain 
the J-averaged FP equation for the energy probability 
density n (E), 


dn{E) 

dt 


2^EP [ DEEn ( E )} 

~LSj (E , Jic) , 


d_ 

dE 


[D E n(E)\ 


( 9 ) 


with effective velocities v E = S E /n (E, J), and vj = 
Sj/n (E, J). The two-dimensional flow in phase-space v = 
(vj,v E ) can be visualized by the streamlines]^] vj/dJ = 
v E /dE, which are derived below from the steady-state so¬ 
lution of the FP equation ( |28| . 

Using the streamlines, we show in Figures |2]|3] the effects 
of the various physical mechanisms. The probability cur¬ 
rent densities are determined by the distribution function (DF) 
of the background stars through the DCs and by DF of the 
test stars. We begin by assuming that a relaxed cusp will 


be appro ximately isotropic (i.e., / (E,J) oc E 1 ^) (Bahcall 
|& Wolf||l976| BW76). The existence of a loss-cone intro¬ 
duces a logarithmic correction, so that the DF is of the form 
f(E,J) oc E 1 / 4 log ( J/Jic) / log ( Jc/Jic) (see Eq. ( |28| ), 
Figures |6] [7] and Hopman & Alexand er|2005| ). Since the DCs 
are not strongly affected by this small anisoptropy, it is con¬ 
venient to assume that the DF of the background is isotropic. 
Therefore, in the calculation of the streamlines we use a DF 
of the form / (E, J) ex /J l/M for the background and a DF of 
the form / (E, J) cx E 1/4 log (J/ J ic ) / log ( J c /Ji c ) for the 
test stars. 

The flow at a point in phase-space is considered j- 
dominated when the streamlines are approximately horizontal 
(i.e., \j \ /j |d| /a) and a-dominated when the streamlines 
are approximately vertical (i.e., j /j <C |d| /a). As shown in 
Figures [2]-|4] the flow is j-dominated, apart for two restricted 


2 This generalizes the simpler situation where stars are only destroyed once 
they reach some high energy threshold, where the loss is expressed instead 
by a boundary condition (cfjBahcall & Wolf|1976f. 

3 The streamlines are immutable under coordinates transformation and 
therefore do not depend on the specific choice of coordinate system. 
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Fig. 2. — Streamlines of the phase-space flow v = (j/j, a/a). All dy¬ 
namical effects apart from resonant relaxation, are included (i.e., two-body 
relaxation, GW emission, GR precession and mass precession). The DCs 
are calculated for a Milky Way-like model (isotropic cusp / ( E) oc E 1 ' A 
with a MBH of 4 X 10 6 ATq, a mass-ratio of M./M* = 5 X 10 5 , and 
total stellar mass M* (rf,) = 2AT. where ry, = 2pc). In addition, 
the probability current densities are calculated assuming a DF f (E, J) oc 
E 1 / 4 (2 J/ J/j) log ( J/ Ji c ) l log ( J r /.J[ r ). which is the steady-state solu¬ 
tion in the presence of a loss-cone. The color map describes the magnitude 
of the DCs. The solid cyan line is the GW separatrix. The dashed cyan line 
indicates the critical sma for EMRIs. 

regions in phase-space (J -A J c and the GW-dominated re¬ 
gion). Therefore, the full flow field can be separated into two 
one-dimensional flows, a fact that will be used below to sim¬ 
plify the analytic treatment. Since RR drives stars only in the 
ji-direction, this separation is enhanced in the phase-space re¬ 
gion where D RR (a,j) > D^ R (a, j), and RR governs the 
dynamics (see Figures |4]-[3| >. 

As shown in Figure [3] in the absence of GR precession 
(i.e., mass precession-only) interior to some sma, RR domi¬ 
nates the dynamics all the way to the loss-cone. In this hy¬ 
pothetical case, the loss rate could be high enough to actually 
empty the cusp close to the MBH. This would then invali¬ 
date the assumption of a single power-low cusp. However, in 
realty, GR precession does play a crucial role in determining 
the dynamics and the steady-state. In fact, due to GR pre¬ 
cession, RR is totally quenched by the adiabatic invariance 
(AI) of the angular momentum. This happens for stars with 
angular momenta below the locus where the precession fre¬ 
quency uip(a,j) = | u>gr + cum | equals the coherence fre¬ 
quency 2tt/T c , where c ogr and cum are the GR and mass pre¬ 
cession frequencies (see AI curve in Figure |4j. Only above 
the AI line can RR be effective. Figure [4] shows a closed con¬ 
tour, somewhat above the AI line, whereRR is faster than NR 
and ther efore dominates the dynamics. As we show in Sec¬ 
tion |3.6| the fact the RR is not effective near the loss-lines 
means that RR does not play an important role in setting the 
steady-state and the loss-rates. 

3.3. Steady state distribution and loss-cone flux 

We assume that the system relaxes in J much faster than 
it relaxes in E. This assumption can be justified by noting 
that J £ [Jic, J C (E )] is bound, whereas E is unbound. We 
therefore assume that the relaxation process is separable: on 
short timescales stars exchange only angular momentum but 


Fig. 3.— Streamlines of the phase-space flow v = (j/j, a/a). All dy¬ 
namical effects apart from GR precession are included (same cusp model as 
in Figure|2j. In this case RR dominates the dynamics all the way to the loss- 
cone. Note that this scenario lea ds to strong depletion of the cusp near the 
MBH (see for example Figurc[T7lBottom). This means that in that region our 
assumption of a power-law steady-state cusp does not hold. The solid black 
line marks the region where RR is effective (stronger than NR). The color 
map describes the strength of the DCs. 


4.0 


3.6 


3.2 


2.8 


I 2.4 


I 2.0 

I 1.6 


1.2 


-2.5 -2.0 -1.5 -1.0 -0.5 0.0 

loglO (J/Jc) 

Fig. 4.— Streamlines of the phase-space flow v = (J /j, a/a). All dynam¬ 
ical effects included (same cusp model as in Figure [2). The solid black line 
marks the region where RR is effective (stronger than NR). The gray line is 
the locus beyond which RR is totally ineffective due to adiabatic invariance. 
The color map describes the strength of the DCs. The solid cyan line is the 
GW separatrix. The dashed cyan line indicates the critical sma for EMRIs. 



not energy and reach their steady-state J-distribution at fixed 
E , and only on a much longer timescale do they reach global 
steady-state in E. This assumption of local equilibrium (i.e., 
in each energy bin the J-distribution is relaxed) is further 
supported by the pattern of the probability current densities 
(Eqs. (|2J, (|3]i) described by the streamlines shown in Figure[5] 
for the NR-only case. The inclusion of RR only strength¬ 
ens this separability. This demonstrates that the motion in the 
.E-direction (a-direction) occurs only at j —> 1, whereas it 
is almost entirely in the j direction at j < 1. The validity of 
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3.6 

3.2 

2.8 

| 2.4 

2.0 

1.6 

1.2 

0.8 


loglQ (J/Jc) 


Fig. 5.— Streamlines of the phase-space flow v = (j / J ■ d/a). Resonant 
relaxation and GW emissions are not included (same cusp model as in Fig¬ 
ure}^}. The color map describes the strength of the DCs. 

this assumption is verified by the excellent match between our 
analytic predictions and the result of MC simulations, which 
do not assume separability a-priori, as shown in Figures [6[- 
[8] In this section we use this separability assumption and 
the fluctuation-dissipation relation (AppendixlE} to derive the 
steady-sate state (E,J) distribution and the flux loss-cone flux. 


In the limit where there is no energy exchange between 


stars, the FP equation can written as (Bar-Or & Alexander 
|20l4) 


dn(E,J,t) 1 d 


d 


dt 2dJ\ JDjj(yE,J " ) dJ 
dSj(E,J,t ) 
dJ 


-n{E,J,t) 


(13) 


which generally follows from the maximum entropy principle 
and in fact provides a necessary test for the validity of the 
DCs (Appendix |E]i. In the absence of a loss-cone (i.e., Jj c —► 
0), the steady-state probability current density Sj (E, J , t) is 
zero, and the local equilibrium distribution is isotropic 


Once the system achieves local equilibrium in J at any E 
(Eq. (17)), the subsequent steady-state in E is obtained by 
solving Eq. 0 for n(E), given Sj(E), 

^ = -5 J(E ) = A {1^ [5ee „ (£)] _ d e „( £ )} , 

(18) 

where R p (E ) = — fj^ Se {E, J) dJ, is the cumulative num¬ 
ber of stars lost through the loss-cone per unit time. Note that 
in steady-state, the probability current density in the J direc¬ 
tion equals the probability current density gradient in E (from 
continuity considerations: the density carried by the J-current 
at fixed E and lost through j R . is balanced by the /(-gradient 
of the total /(-current). 


3.4. Steady state distribution for two body relaxation 

We now show that in the case of two-body relaxation, the 
solution of the energy FP equation (Eq. ( [18) ) with non-zero 
flux can be approximated analytically to derive the steady- 
state density distribution and the plunge rate. Since the plunge 
rate is small compared to the relaxation rate, the energy distri¬ 
bution asymptotes to the zero-flux (i.e., no plunges) |Bahcall| 
& Wolf (j 1976}) power-law solution. 


For two-body relaxation, d p R ( E,j ) asymptotes to a finite 


33 

value d° NR (E) = df R (E, j = 0) as 


the contribution to the current density, 


5 j -A (R Si 
ity, Sj[E), 


Since most of 
reflects the 


value of djj at small j (Eq. ©), it can approximated by 
d%{E), so that 


and 


Sj(E) 


~n{E) 

~n(E) 


_ 2d% R (E) _ 

log (J c 2 /Jf c ) ~ 1 + JfJJ 2 

log (J c /Ji c ) ’ 


n {E, J) 


n (4 

44 


2 J log ( J 2 /Jf c ) 

J c log UcM 2 c) - 1 + 3lc 
2 J log ( J 2 / Jf c ) 

J 2 c log (J c 2 /4) ' 


(19) 


( 20 ) 


n (E, J ) = 2 J/J 2 c (E)n{E). (14) 


For a finite J/ c , it follows from the separability assumption, 
that the probability current density is non-zero and indepen¬ 
dent of J. Therefore, from Eq. ([13) we obtain 


JDjj(E,J ) 


d_ 


1 

J 


n(E , J, t) 


-2 Sj (E) . (15) 


Since dE R scales as Dee/E 2 = T E X (Appendix |Ci, it is 
convenient to represent Sj(E) explicitly in terms of the en¬ 
ergy relaxation time, T R , as Sj(E) = —n (E) \{E) /Te{E), 
where 


x(E) = d° NR (E)/ log (J c /J lc ) = 2d% R {E)/\og{E lc /E) , 

( 21 ) 


By integrating over J and using the normalization n (E) = 
n(E, J)dJ, we obtain 


and 


Sj (E) 


~n{E) 


If 1 1 ~ j 2 dj 

/ 4c djj(E.j) j 


(16) 


■ (E ’ J) — SAE) Lfwd 3 ' <17) 


were djj (. E,j ) = Djj (E, J ) /J 2 . 


expresses the logarithmic suppression of the flux due to the 
decreasing size of the loss-cone away from the MBH, and 
where E) c = GM, /32r g corresponds to the limit J c = J; c 
for J[ c = 4r g c and for Keplerian energy. Note that this is 
not the true innermost stable circular orbit, but rather a for¬ 
mal extrapolation of the approximations adopted here, used 
for normalization only. Here we are interested in stars with 
E Ei c where 1/ log (Ei c /E) is small. The cumulative 


Since dff 1 = 2 jdf R for j —> 0 (e.g., Shapiro & Marchant 1978 


7 7 J J v ■f’ l 

AppendixjcJ and 2jdf rt = djjd^- R (AppendixjE). 
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plunge rate Rp = — f SjdE can then be approximated as 


R 


p 


\n(e) 


X(E) 

T e 



X(E') 

T e (. E') 


N(E') 


dE’ 

~W 


( 22 ) 


where N (E) is the number of stars with energy larger than 

E. 

For an infinite isotropic cusp n (E) oc E p ~ 5 / 2 where 0 < 
p < 1/2 (to ensure the DCs are finite), the J-averaged DCs 
are ( |Bar-Or et al.|20l5) 


D e = 


Dee — 


4(1 — Ap) log A N(E) 
1+p Q 2 P{E) 
16 (3 - 2p) log A N(E) 
(l+p)(l-2p) Q 2 P{E) 


E(xE p+1 , (23) 
E 2 (xE p+2 , (24) 


and the plunge rate is 



R 


p 


3-2 p 
3 — 4 p 


N(E) 


X(E) 

Te 


(25) 


Since \ is nearly constant in E for E << Ep., it can be approx¬ 
imated in that limit by evaluating it at E m i n = GM,/ 2a max . 
In that case, Eq. © has an analytical solution. 


(4 p - 3) (1 - 4 p) 
4(3-2 p) 


(26) 


which connects the current to the power-law exponent of the 
cusp. The physical branch of the solution i^j 

P=\ (-Vx 2 +8* + 1 - X + 2) « j (1 - 5x) ■ (27) 


Since for p = 1/4, d, 0 j:) « 1/(8 Te) (Appendix |c 1 , it fol- 


lows that y ~ _ 1/ [8 log (1 /ji c )\ <C 1, and so the Bah 


call & Wolf (1976 1 solution p = 1/4 (i.e., y —» 0) is a 
reasonable-approximation in the E —> E m - ln limit, where 
ji c << 1. Thus in steady-state the energy (or sma) distribu¬ 


tion is n (E) oc E 9 / 4 (or n (a) oc a 1 / 4 ). Using Eq. (|20j we 
obtain the (E, J) steady-state distribution 


Fig. 6.— The energy distribution as function of energy shows a good agree¬ 
ment between the analytic BW76 cusp solution and MC re sults . The critical 
sma, clgw(Egw) as defined by the sparatrix (see Section [33) is shown by 
a dashed line. 



/Z? T\ ^ ]\J ( Z? ^ ( J? / J? 9 / 4 

' ( ’ J) ~ 4 ( min) Jl log ( J c / J lc ) ( 1 min) 


N(a) 


and a steady-state plunge rate 

p . 1_i_ 

32 log {1/jic (a)) T e (a) 

where the energy relaxation time is ( |Bar-Or et al.|2013] > 


Te (a) = —Q 


P ( a ) 


64 N (a) log Q 


(28) 


(29) 


(30) 


In the limit E —>• Ei c , x(E) can no longer be approxi¬ 
mated as fixed nor as small, and the power law solution breaks 
down. We now argue that this power-law approximation is 
valid for almost the entire range of MBH masses and their 


5 This generalizes the analytic BW76 solution (p = 1/4,\ = 0), which 
applies for a power-law DF in steady-state with a constant .E-current (which 
then must be zero). Here, the “leakage” of stars through the loss-cone at all E 
(X > 0) implies a non-constant current, which allows flatter cusp solutions 
withp <1/4. 


Fig. 7.— The angular momentum distribution in a bin centered around 
a = 0.4a max . T he M C data shows the expected logarithmically-suppressed 
distribution (Eq. |28j) compared with the isotropic one. 


host galactic nuclei. Define apw as the minimal sma where 
X and 31ogx/<91og.E are still small enough for this approx¬ 
imation to hold. Since 91ogx/<9k>g.E oc x oc 1/log (1 /ji c ) 
where ji c oc yjr g /a , it follows that clew <x M.. We are 
interested to resolve the dynamics close to the MBH in order 
to obtain a relia ble e stimate for the EMRI event rate. As we 
show in Section 13.51 the rate is determined by the dynamics 
near the critical GW sma, acw ■ Thus it is sufficient do show 
that clew < (I cw^ since then power-law approximation is 
valid over the entire range of relevant radii. The MC results 
shown in Figures [6] and [8 for the case of M. = 4 x 1O 6 M 0 , 
demonstrate that our analytic estimates (Eqs. @, @) for 
the steady-state energy distribution and plunge rate reproduce 
the simulated results at least down to acw- This holds also 
for more massive MBHs. The M/a relation and the scalings 

clew oc M. and acw oc (log Q) -4 ^ 5 (Eq. ( |A1Q| )) imply 

clew/ci-GW oc (logQ) 4/5 M. 2/,J and therefore clew «Gtv 
up to M, ~ 10 9 Mq. 



































10 


Bar-Or & Alexander 



Fig. 8.— The plunge rate as a function of semi-maj or a xis. The MC data 
(without RR) agree with the analytic expression (Eq. (46p . The MC results 
demonstrate that the contribution of RR to the plunge rate is small. The 
critical sma, ciqw (vertical dashed line) as defined by the separatrix fsce |3.4| 
agrees very well with the sma (vertical solid line) where the inspiral rate, 
R\ ot (dash-doted lines), equals the plunge rate R p ° Gw in the absence of 
GW emission. 


3.5. EMRI event rates 

So far we ignored the contribution of GW emission to the 
dynamics. Compact objects can withstand the tidal field of 
the MBH. When on eccentric orbits, their orbital decay by 
the emission of GWs can be faster than the diffusion of angu¬ 
lar momentum due to the stochastic perturbations of the stellar 
background. In that case, they inspiral gradually all the way 
dow n to the innermost stable circular orb i t (ISCO) a s EM- 
RIs ( [Peters & Mathews|1963[|Peters|1964l|Gair et al.|2006) , 

instead of plunging directly into the MBH with J < (Fig- 
ure[T]i. The GW signature of plunges and inspirals is very dif¬ 
ferent. The low mass of the compact objects generates weak 
signals, well below the noise. Plunges result in short, very 
hard to detect broad spectrum GW flares. In contrast, EMRIs 
are of special interest since their quasi-periodic signal can be 
integrated and detected against the noise if the waveform is 
approximately known. 

In the absence of GWs (Figure [5]), the streamlines are ap¬ 
proximately constant in a. In contrast, GW emission diverts 
the streamlines to tracks that are almost parallel to the loss- 
cone (i.e., nearly constant J) in the phase-space region where 
GW dominates the dynamics (Figure[2j. The outermost inspi- 
raling streamline separates phase-space into two distinct re¬ 
gions. Above this separatrix all streamlines are plunges, while 
below it all streamlines are inspirals (Figure [4]). The conti¬ 
nuity equation (Eq. ( fT2] i) implies that the probability current 
in steady-state is constant along a streamline bundle. Since 
the streamlines in the GW-dominated region below the sep¬ 
aratrix originate in phase-space regions where the density is 
much higher, the small depletion due to EMRI losses is not 
expected to affect the density at the origin of the streamlines. 
The EMRI rate can therefore be estimated by identifying the 
terminal point (a p , ji c ) of the plunge streamline (without GW) 
corresponding to the separatrix. This is the effective critical 
sma for EMRIs, clow- The EMRI rate is then obtained by 
integrating the differential plunge rate in the absence of GW 



Fig. 9.— The plunge ratio R GW /Rp° GW , where R GW and R p oGW 
are the cumulative plunge rates obtained with and without GW dissipation. 
The crit ical sma, acw (vertical dashed line) as defined by the separatrix 
(see |3.5fr agrees very well with the sma (vertical solid line) where the inspiral 
rate, R\ ot , equals the plunge rate R Gw in the absence of GW emission. 

emission from a isco to clcw* 

f E isco 

Rl ot = - / Sj {E) dE = R p ( a GW ) ■ (31) 

J e gw 

Thus, for a BW76 cusp the EMRI rate is 

R tot = j>_1_ N (acw) 

32 log {1/jic {acw)) Te (o-gw) 

_ acw log iX/jic (a max)) u noGW („ \ /qox 

“ a max log (1 /j lc (a GW )) p [ j 

which is approximately linear in acw- The value of acw 
is determined by solving the streamline equation dE/dJ = 
Se/Sj with the boundary condition that the streamline tra¬ 
jectory reaches at E max . The probability current densities 
are estimated by assuming that Sj is constant in J and Se 
results only from GW dissipation. This means that in the ab¬ 
sence of GWs, that streamline is constant in E and acw can 
be estimated by taking the value of a at J J; c . The exact 
value of acw depends on the GW emission approximation 
used and is calculated in Appendix [A] As shown in Figure [8] 
this definition of acw is indeed a good approximation to the 
sma where the plunge rate (without GW) is equal to the inspi¬ 
ral rate and can be used to predict the inspiral rate in the MC 
simulations (Figure [l7]i. As expected, inside a gw the plunge 
rate with GW decreases relative to the plunge rate without 
GW (See Figure 0. 


3.6. Effect of Resonant Relaxation 

Due to the long coherence time, resonant relaxation is a 
much more effective process than two-body relaxation. How¬ 
ever, in regions of phase-space where in-plane GR precession 
is faster than the coherence time, j becomes an adiabatic in¬ 
variant and the RR process is quenched. RR it is therefore lim¬ 
ited to a small region of phase-space (see Figureffli. The l ocus 
where in-plane precession quenches RR by AI (Section 14.1 1 
defines the outer envelope of the region where RR may be ef- 
ficient relative to NR. The region where RR dominates the 
dynamics even on long timescales is where by the ratio of the 
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2nd order DC^jexceeds unity, i.e., Djj RR /Djj >RR > 1 (See 
Figure [4 1 . The typical phase-space configuration is shown in 
Figure ]T the region where RR dominates is detached from 
the loss-lines; NR is required for the stars to evolve towards 
them, and therefore the slow NR timescale remains the bot¬ 
tleneck for the loss-rates which are mostly unaffected by RR 
(see Figure [TT}. This can be shown formally by re-estimating 
the probability current density in the presence of RR. 

For a smooth noise, the RR diffusion coefficient can be 
written as ( |Bar-Or & Alexander|2014j > 

= 2 T coh iyje~ 4 ^o^\ (33) 


where vj = \/ ( t 2 t ) is the residual torque in the J direction 
(see Appendix ] d]i 

vj » (a/2)J c v a /Q, (34) 

and j 0 (a) is the AI locus where the GR precession frequency 
i/Gr (a, j) equals to the coherence time 


Jo 


= \/2nT c UGR ( a,j = 1). 


Here we adopt 


T c = V^]2v~ l (2a, x/ITs) 


where v p is the combined mass and GR precession. 
The combined diffusion coefficients are 

Dj = D^ r + D^ r , 


D JJ = Dy j R + D' JJ , 


■)RR 


and using Eq. ([T6l>, the flux is given by 


Sj (E) = —n (E) / 


1 ~j 


dj 


Djf/J/ + D RR /J/ j 


(35) 


(36) 


(37) 

(38) 


(39) 


Since D R R is rises up to some maximal value before sharply 
drops as it approaches j —> jo, we can approximate the RR 
DC as 


rjRR 


d RR^-j)e 47rj ° j>j„ 
0 j < j n 


(40) 


where d° RR = 2T co y i i'j ( j = 0) and the maximum 
of d RR (E,j), occurs at j m , given by j^/jg = 
167T (1 — j rn ) /jo ~ 167T (1 — jo) /jo- The differential flux 
is therefore given by 


Sj R (E) 

W) 


1 ~ j 2 dj f 1 l~j 2 dj 

4 = d ?J R ■> Jj m d^ R + df R j 

1 _ Xrr log (l67r (1 — j 0 ) jp) 

1 + Xrr 5 log (ji c ) 


-1 


(41) 


where \RR = d° RR e-' l ' K: i° /d° NR - As shown in Figure[lo] this 
analytic approximation reproduces the MC results. 

The small effect of RR on the loss-rates can be estimated 
by integrating Eq. ( |4Tj ) over the relevant region 

P -E'max 

R rr {E) = - Sf R {E') dE '. (42) 

Je 


6 The transformation of the DC from J to j = J/J c is = 

Djj,nr/J l + U 2 /4)D E e/E 2 +jD EJ /J c E (Appendix[cj. 



Fig. 10.— The ratio between the pr oba bility current density in J with 
and without RR, as function of a (Eq. ( |4 1 ^ ), which expresses the differen¬ 
tial plunge rate. The limited increase in the ratio over its asymptotic value of 
1 in the a —¥ 0 and a —>■ oo limits expresses the fact that the contribution of 
RR to the plunge rate is small. 


4. MONTE CARLO MODELS 

We complement and validate our analytic study of the rel¬ 
ativistic loss-cone by numerically evolving the FP equation 
in both E and J using a MC procedure (described in Ap¬ 
pendix [B]). Unlike the analytic treatment, this procedure does 
not assume that the evolution in J can be decoupled from that 
in E. The advantages of the MC method over direct A'-body 
simulations are the high degree of flexibility it offers for iso¬ 
lating and studying the different mechanisms that affect the 
dynamics of the loss-cone, the ease of including additional 
physical effects and of modifying the initial and boundary 
conditions, and importantly, its scalability to systems with 
a realistically large number of stars. In Section [5] below we 
employ the MC procedure to calculate the phase-space den¬ 
sity and rates of relaxed galactic nuclei, and in particular, that 
of a Milky Way-like nucleus with a 4 x 10 6 M ( . } MBH and 
N ~ O(10 7 ) stars on its radius of influence, ~ 2 pc. Such 
nuclei are considered archetypal for future space-borne mis¬ 
sions to detect low-frequency gravitational waves from inspi- 
raling compact objects. 

We begin here by validating the MC procedure. We study 
the dynamics in the restricted case where E remains fixed, 
which allows a direct test of the impact of a diabatic invari¬ 
ance on the long-term dynamics (Section 14.11. We then com¬ 
pare rate results from our MC procedure in both E and J with 
the currently available results from direct post-Newtonian 
IV- body simulations of small-TV (N < 100) systems (Sec- 
tion |4.2| ). 


4.1. j-only Monte Carlo simulations 

The maximal entropy limit (Appendix [E]) provides a basic 
test for the physical validity of the DCs and of the MC pro¬ 
cedure for evolving the Fokker-Planck equation. The proba¬ 
bility density of a closed system with zero net angular mo¬ 
mentum must asymptote to the maximal entropy solution 
dN/dj = 2j. Experimentation shows that this is a sensi¬ 
tive test of both the functional form of the DCs and the de¬ 
tails of the MC procedure, in particular the implementation 
of the boundar y cond itions. We verify the maximal entropy 
limit in Section [4. 1.1 1 In the absence of NR (for example on 
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Fig. 11.— The convergence of a MC y-only NR diffusion simulation to 

the maximal entropy solution dAf/dJogi n 1 = 2 In 1 Of 2 . . __ 

timescales <C 1 nr). a relativistic system that is subject to RR 

with a smooth background noise should display adiabatic in¬ 
variance (AI) in the form of a sharp drop in the phase-space 
density below some small value of j where the GR preces- 
sion perio d falls below the coherence time (Bar -Or & Alexan-1 
|der|2014| l. RR with non-smooth background noise is not ex¬ 
pected to display such an AI barrier (the “Schwarzschild bar- 
Merritt et al.|2011 1 . We demonstrate that o ur MC pro¬ 


cedure reproduces this behavior in Section |47L2| Finally, we 
study the rea listic c ase where NR smears the RR-generated 
AI in Section |4.1.3| and also show how this smearing appears 
in the unrestricted case where both a and j evolve. 

Since the j —► 0 limit is of special interest, it is efficient 
to use logarithmic bins to collect statistics on the phase-space 
density. In that case, it is more useful to represent the den¬ 
sity as dN/dlog jj^vhere the maximal entropy solution is 
dN/ d log j = 2 j 2 (or dN/d log 10 j = 2 log(10)j 2 ). 


4.1.1. NR only 

Figure 11 shows the j-PDF at a = 


flmax/4 for ana = 
= 10 4 a min , after time t = lOOTg 
). Near-c omple te convergence is al- 


7/4 cusp with a max 
(Te = [E 2 /A(E 2 )] a 
ready reached at T > Te, (Section|4. 1 .3k. The convergence to 


the expected maximal entropy solution is apparent, although a 
bias toward a somewhat steeper slope for j < 0.1 is observed. 
In addition to the ct = 7/4 case shown in Figure ITT] simula¬ 
tions with other values of a or of a max /ctmin confirm that the 
maximal entropy solution holds generally. 

4.1.2. RRonly 

Figure [12] shows the j-PDF for the RR-only case. The 
MC code reproduces the AI barrier for_smooth noise at 
jo = V T c“Gr/2 it (Bar-Or & Alexander 2014;), while non¬ 
smooth noise asymptotes as expected to the maximum en¬ 
tropy solution—a demonstration that this limit does not de¬ 
pend on the nature of the relaxation process. 

4.1.3. NR+RR 

The presence of NR erases the AI cutoff in the j-PDF on 
timescales approaching or exceeding the NR timescale (quan¬ 
tified by the energy diffusion timescale Te). This is demon¬ 
strated in Figure [T3] which shows a sequence of j-only MC 

7 The decreasing size of the bins in linear space leads to a misleading 
graphical representation of dNjdj when the statistics are low, as is the typ¬ 
ical case at low-J, since the normalized bin density AN/Aj diverges for 
AN > 1 as Aj -s- 0. 



Fig. 12.— Adiabatic invariance and the convergence of j RR to the max¬ 
imal entropy solution dlV/d log 10 y — 21nl0y 2 for two different back¬ 
ground noise models: non-smooth noise (C°) with an exponential ACF (top) 
and smooth noise (C°°) with a Gaussian ACF (bottom). The predicted char¬ 
acteristic position of the AI front at jo is marked by the vertical line. An 
approximate upper limit on the density (1 count per bin) was estimated for 
empty bins (triangles). 

simulations that include both NR and RR. All the simulation 
runs had a fixed duration T s i m = 100 Te, which kept the num¬ 
ber of binned points, and hence the statistical sampling fluctu¬ 
ations, fixed (the MC values are sampled every At = 1 T e ). 
However, the effective NR timescale was artificially extended 
to T' e = X d T e (X n > 1), so that T sim = (100 /X D )T' E , 
and Xe> was varied from Xp = 10 3 (T s j m = 0.1T^) down to 
X E = 1 (T s ; m = lOO Tp): t he larger X E , the less significant 
is NR over T s ; m . Figure 
^sim = T' e , and T s ; m = 
smeared already when T E 
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shows results for T s ; m = 0.1X^ 
Bfl T' e . The AI cutoff is substantially 
= 0.17'/ (compare Figure [l3| top 
left with Figure [12] right), and the AI remains only as a mod¬ 
erate steepening of the slope below jo for T s j m = T' f . For 
T s im = 1007),-, the j-PDF is almost indistinguishable from 
the case of NR-only. 

This trend is evident also in the general case where both a 
and j are free to evolve, as shown in Figure [14] On timescales 
of order of the RR relaxation time, but mucFshorter than the 
NR timescale, the stellar trajectories are bound by the AI line. 
However, on longer timescales, NR drives stellar diffusion 
across the AI line and beyond. The existence of a persistent 
Schwarzschild Barrier with a locus ass(j) oc j~ 2 /( 5 ~ a \ as 
suggested by MAMW11 (Eq. 35 there), is neither supported 
by our analysis nor observed in our MC simulations. 

4.2. Comparison with N-body simulations 
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-3 - 2.5 -2 - 1.5 -1 - 0.5 

l°8ioi 


Fig. 13.— The suppression of adiabatic invariance by NR and the convergence of j-only MC simulations with NR and RR with a Gaussian ACF noise to the 
maximal entropy solution d N/d log 10 j = 2 In 10j 2 . All MC runs lasted a fixed time T s [ m = IOOTe, were Te is the energy diffusion timescale at a max , but 
the NR relaxation time was artificially extended to T' E = XqTe (Xp > 1) so that T s [ m = (100 /Xe)T' e (the larger Xp, the less significant is NR). The 
predicted characteristic location of the AI barrier at jo is marked by the vertical line. The best fit power-law to the j-PDF slope is also shown (green line). An 
approximate upper limit on the density (1 count per bin) was estimated for empty bins (triangles). Top left: T s i m = O.lTg,. Top right: T s i m = T' E . Bottom left: 
T s im = lOOTg = IOOTe- Bottom right: NR-only MC simulation for comparison. 



Fig. 14.— MC snapshots of the trajectories of individual stars at different times, corresponding to increasing fractions of the energy relaxation timescales Te 
at a max . On short timescales, t < Tg, stars do not cross the AI line (solid line), while on longer timescales, t — >■ Te, NR progressively drives stellar diffusion 
across the AI line to the entire available phase-space. The MC simulation assume a MBH of 10 g Mq and a cusp of 50 stars of 50M® each with initial conditions 
drawn from an isotropic cusp with a. = 7/4 and a ma x = 10 mpc. We also plot for comparison the Schwarzschild barrier for this cusp model (dashed line), as 
suggested by MAMW11. 
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log, 0 (<rn/dlgj dlga) 



- LSO 

- GW 

v insplral 



logioO) 



Fig. 15.— The 2D loss-cone phase-space density in a MAMW11-like cusp 
model with mass-precession coherence time and Gaussian noise, calculated 
by high phase-space resolution MC simulations. Top: Only NR. Middle: Full 
model with all the dynamical processes. Bottom: Same as full model, with 
GR precession switched off. The black contours denote the loci where RR 
dominates NR, D RR /D^ R = 1,10,100. To avoid clutter, only 1% of the 
plunge and inspiral terminal track points (circles and triangles) are displayed. 


low the stellar number to fall with time as stars are lost into 
the MBH. In addition, the MAMW11 and BASH models 
all have an initial n*(r) ex r -2 cusp, which is away from 
the BW76 steady-state configuration of a single mass popula¬ 
tion, n*(r) oc r - '/ 4 . Thus, these A-body simulations always 
remain out of steady-state due to relaxation, expansion and 
stellar loss. The loss-rates of the MC (when its parameters 
are matched to initial state of the A-body simulations) should 
therefore be compared to the initial loss-rates of the A-body 
simulations. To further reduce this incompatibility, we mod¬ 
ified our MC procedure to reproduce the initial conditions of 
the A-body simulations by introducing the test stars into the 
interior of the cusp according to an ?i* ex r -2 probability den¬ 
sity. 

Figure p~5| shows the phase-space density derived from MC 
simulations for the MAMW11 cusp model (see details below), 
which is used here for comparison with A-body results. Re¬ 
sults are shown for the case of Newtonian dynamics without 
RR, and for the case where all dynamical effects are switched 
on (Newtonian dynamics, GR, NR and RR). Also shown are 
the endpoints in phase-space of a representative fraction of 
the plunge and inspiral events. Figure [15] and Table [2] demon¬ 
strate that GR precession plays a critical role in making EM- 
RIs possible; in its absence, RR remains unquenched at all 
low values of j > ji c and therefore stars are rapidly driven 
to plunging orbits before they can reach the EMRI line and 
inspiral by the emission of GW. However, when GR preces¬ 
sion is included, the region where RR dominates over NR is 
restricted to regions that are far from the loss-lines (black con¬ 
tours, Figure [15] top right). This creates a bottleneck in the 
flow from a max and j ~ 0(1) to the loss-lines, where the 
orbital evolution is driven by slow NR, and therefore the ef¬ 
fect of RR on the loss-rates in steady-state is not large. This 
near-independence of the steady-state loss-rates from RR is 
analyzed in Section [3] The MC results also show that mass 
precession cannot play a similar role, since it becomes effi¬ 
cient only for j —> 1 and a —>■ r ? l orbits. 

Computational costs limited the MAMW 11 simulations to 
a non-realistic cusp of only A* = 50 heavy objects of mass 
M* = 50 M 0 each around a MBH of mass M. = 10 6 Mm, 
extending between a m ; n = 10" 3 mpe to a max = 10 mper] 
The stars were initially set on stable orbits, isotropic in ori¬ 
entation and eccentricity. GR was introduced to the equa¬ 
tions of motion perturbatively, up to post-Newtonian (PN) 
order PN2.5. Orders PN1 and PN2 contribute only to the 
in-plane (Schwarzschild) periapse precession, while order 
PN2.5 contributes only to dissipative GW emission. By selec¬ 
tively switching on or off the various PN terms, the A-body 
simulations tested the cases of Newtonian gravity (all PN 
terms switched off), no GR precession (only the PN2.5 term 
switched on), or full perturbative GR (all PN terms switched 
on). Table [2] compares the loss-rates for the corresponding 
MC and A-body simulations, as well as for an artificial model 
that can only be realized by the MC method, a Newtonian case 
where RR is switched off. The MAMW 11 loss-rates were re¬ 
ported as function of time, so it is possible to extrapolate to 
t —> 0 and obtain lower or upper limits on the rates. The 


Matching results of MC simulations to results from direct 
A-body simulations is not straightforward. The MC proce¬ 
dure enforces boundary conditions at a max and assumes an 
approximate steady-state background, whereas the A-body 
simulations of MAMW 11 and BASH provide a max only as 
initial conditions (the cluster subsequently expands), and al- 


8 These constraints lead to atypical dynamical properties in this model. 
(1) Because TV* is so small that y/ TV* ~ 0(TV*), the difference between 
the mass precession coherence time and the self-quenching coherence time 
is small, and the two are hard to discriminate. (2) Since the NR relaxation 
time is T^r ~ Q 2 P/TV* log Q. while the mass-precession RR timescale 
is Trr ~ QP, Trr/Tnr <K 1 everywhere in the cusp, that is, RR is 
atypically efficient. 
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BAS 14 rates where reported only in the average. 

Table|2]shows that the MC loss-rates for the full GR models 
are quite similar, irrespective of the RR noise model, and that 
they are also similar to the rates predicted in the artificial case 
where RR is switched off. The MC loss-rates are somewhat 
higher than those derived from the iV-body simulations. The 
MC model with smooth (Gaussian) noise provides a better fit 
to the A'-body results, while the coherence models are vir¬ 
tually indistinguishable, with only a marginal preference for 
mass precession. 

Overall the MC loss-rates are in agreement with the N- 
body simulations; the systematic trends can be explained by 
the differences between the computational techniques and the 
physical assumptions. Compared to the MAMW11 TV-body 
simulations, the MC plunge rate is 3.4 higher and the in¬ 
spiral rate is 1.7 higher; compared to the BAS14 IV-body 
simulations, the MC plunge rate is 1.5 ± 0.3 higher and the 
inspiral rate is 1.8 +o'q higher. Since the incompatibility of 
the MC and IV-body treatment of the boundary condition at 
a max results in a lower stellar density in the A'-bodv simu¬ 
lations, the fact that the A'-body loss-rates are systematically 
lower is to be expected. The same systematic trend is also 
seen in models where GR is switched off (i.e., no GR pre¬ 
cession, no GW). When only GW is switched on but GR 
precession is switched off (i.e., no GR quenching), the MC 
inspiral rate is zero in agreement with MAMW11. An addi¬ 
tional difference between these studies is the plunge criterion. 
MAMW11 used r < 8 r g , BAS 14 r < 6 r g and here the cri¬ 
terion was based directly on the angular momentum, J < Ji c 
(where J was evaluated in the Keplerian limit). As noted 
by Gair et alT] (|2006|>, plunging orbits (i.e., parabolic orbits 
with J = Jic = At!M/c) correspond to Keplerian orbits with 
periapse r/ c = 8 GM/c 2 , or to relativistic orbits with peri- 
apse Tj c = AGM/c 2 . This can explain some of the systematic 
differences in the rates, since if r; c over-estimates the true 
value, stars that should inspiral plunge prematurely, thereby 
biasing the rates to too-high plunge rates and too-low inspiral 
rate. Conversely, if r; c under-estimates the true value, a too- 
high inspiral rate will follow. We believe that this explains the 
discrepant non-zero inspiral rate that BAS 14 find for the case 
where GR precession is switched off (i.e., Keplerian dynamics 
where ri c —> 8 r g is the correct limit). Our approximate an¬ 
gular momentum plunge criterion for parabolic orbit applies 
generally in both the Keplerian and relativistic regimes ( |Gair 
et al. 120061 , unlike the periapse criteria used by the A'-body 
simulations. This may explain why the MC inspiral rates are 
somewhat higher than the A'-body results, which assume a 
too-high value of r; c for the relativistic regime. 


5. LOSS RATES 

5.1. The Galactic Center test case 

The MBH in the center of the Milky Way and the stars 
and compact objects around it are a system of particular rele¬ 
vance, both because it is uniquely accessible to observations, 
and can therefore place constraints on dynamical models and 
theories Alexander (2005J), and because planned space-borne 
GW detectors with 0(10° km) baseline will be optimally sen¬ 
sitive to GWs emitted by a mass orbiting a 10 6 — 10' Mq 
MBH near the last stable orbit (Amaro-Seoane et al. 


2007). 


Therefore, although it remains an open question whether the 
Galactic center (GC) is surrounded by a high density relaxed 
cusp of stellar remnants (see review by |Alexander|201 1) , and 
although the rate of GW events from the GCitself is expected 


TABLE 2 

The plunge and inspiral rates in MAMW 11 -like cusp models 


Method 1 

Processes 2 

T c 3 

Noise 4 

Plunge 5 

Inspiral 5 

MC 

All 

SQ 

E 

0.96 

1.8 

MC 

All 

SQ 

G 

0.71 

1.5 

MC 

All 

M 

E 

0.92 

1.8 

MC 

All 

M 

G 

0.71 

1.4 

MC 

No RR 

— 

— 

0.68 

1.5 

MC 

No GR 

M 

G 

110 

— 

MC 

No GR prec. 

M 

G 

no 

0 

MC 

No mass prec. 

M 

G 

0.71 

1.4 

NB1 

With GR (t -A 0) 



0.2 

~ 0.9 

NB1 

No GR it -A 0) 



> 20 

— 

NB1 

No GR prec. (t —»■ 0) 


> 20 

< i 

NB2 

With GR 



0.5 ±0.1 0.8 ±0.1 

NB2 

No GR 



28 ± 2 

— 

NB2 

No Gr prec. 



26 ± 2 

4.3 ±0.6 

1 Method: MC = Monte Carlo, NB = iV-body (1: MAMW11, 2: 

BAS 14). 

2 Processes: 

All includes NR, RR, GW 

Gair et al. 

2006 

, mass prec., GR prec. 

6 Coherence time: M = Mass prec., SQ 

= Self-quenching. 


4 Noise model: G = Gaussian noise, E = 

Exponential noise. 


5 Event rates in units of 10 - 6 yr -1 . 






to be small (but see Freita g|2003 i, the Galactic MBH SgrA* 
represents the archetypal cosmic GW target. 

We adopt here a simplified model of the GC consisting of 
an MBH mass of M, = 4 x 10 6 Mq, surrounded by a steady- 
state a = 7/4 cusp of equal mass stars of either M* = 1 Mq 
or 10 Mq, extending between a m j n = 2 x 10 -4 pc (ai n = 
0.1a min , see Appendix |Bji to a ou t = a ma x = = 2 pc with a 

total stellar mass of M*7V*(< rh) = 2M. inside the radius of 
influence. Test stars are injected into the nucleus with initial 
sma ao = a out , and with isotropic j 0 . 

Figure fl6| shows the steady-state configuration and loss- 
rates for aGC model with M. k = 10 Mq and a smooth back¬ 
ground noise whose coherence is set by mass precession. As 
expected from the fact that the RR-dominated region in phase- 
space is well separated from the loss-lines, the steady-state 
phase-space density is very close to the simple NR-only solu¬ 
tion. 

Table 


16 explores the implications of varying some of the 


assumedprocesses for the loss-rates: the mass of the cusp 
stars (1 Mq or 10 Mq), the nature of relaxation (NR only, or 
NR and RR) the noise model (White, exponential, Gaussian), 
the coherence mechanism (mass precession, self-quenching), 
or the GW dissipation approximation. 

The uncorrelated (white) noise model for the resonant 
torques, which is equivalent to the assumption that ugr —> 0, 
results in very high plunge rates as strong RR rapidly drains 
the cusp, and as a result the EMRI rate drops to zero. In con¬ 
trast, for all other RR models, irrespective of the assumptions 
about the nature of the noise or the coherence mechanism, GR 
precession suppresses RR to the extent that rates are very sim¬ 
ilar to those derived for the non-physical “NR-only” model: a 
plunge rate of T p ~ (6 — 9) x 10~ 4 yr _1 , and an inspiral rate 
of Tj ~ (1 — 3) x 10 _6 yr. We find that the more sophisti¬ 
cated GW dissipation estimate of|Gair et al.|(j2006j) results in 
inspiral ra tes tha t a re a factor ~ 2 higher than the estimates 
by Petersl (| 1964j) or Hopman & Alexander ([2006a]). We con- 
clude that to within a factor of ~ 2, our rate estimates for 
relaxed steady-state cusps are robust to variations of the phys¬ 
ical assumptions. 


5.2. Scaling with the MBH mass 
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Fig. 16.— The 2D loss-cone phase-space density in a Milky Way-like cusp 
model with mass-precession coherence time and Gaussian noise, calculated 
by high phase-space resolution MC simulations. Stars / stellar mass BHs of 
10 Mq are assumed. Top: GR precession included. Mass precession limits 
the efficiency of RR beyond ~ 100 mpc, while GR precession limits RR 
below the AI locus (gray line). RR is faster than NR only well away from 
the loss-cone, inside the black contours (equally fast at the outer contour, 10 
times faster at the inner contour). Bottom: When GR precession is artificially 
switched off, RR remain effective all the way down to the loss-cone and is 
faster than NR below a ~ 100 mpc. As a result, stars are driven to plunge 
trajectories well before they can lose enough energy by NR to reach the GW 
loss-line. A central, strongly depleted cavity is formed, and the EMRI rate is 
completely suppressed. 


The MC simulations can be used to validate a simple ana¬ 
lytic model for estimating the loss-rates and their dependence 
on the parameters of the galactic nucleus, which is based on 
identifying critical values of the sma, a c , below which the 


probability of a star to cross the loss-line is 0(1) (Lightman 
|& Shapiro] 1977 Hopman & Alexander 20051 . The loss-rate 
is then i' oc N(< a c )/'lNR(a c ), where the proportionality 
factor includes the suppression of the density near the loss- 
line. For plunge events a c ~ CAfry,) ([Lightman & Shapiro 
1977| >, while for GW inspiral a c ~ a G w fh, the max- 
imum of the GW line (Section [3}. Figure [I] shows that the 
region of phase-space where RR dominates the dynamics is 
well separated from the loss-lines, is well below ry, and well 
above a G w- The timescale relevant for estimating is therefore 
that of NR and not RR. 

In order to estimate the integrated cosmic rates of EMRIs 
or tidal disruption flares, it is necessary to scale the loss-rates 


TABLE 3 

The plunge and inspiral rates in Milky Way-like cusp models 


M* 1 

Processes 2 

T c 3 

Noise 4 

Plunge 5 

Inspiral 5 

1 

NoRR 

— 

— 

730 

3.1 

1 

GW1 

SQ 

w 

16000 

0.0 

1 

GW1 

SQ 

E 

860 

3.3 

1 

GW1 

SQ 

G 

880 

2.3 

1 

GW1 

M 

W 

930 

0.0 

1 

GW1 

M 

E 

840 

3.2 

1 

GW1 

M 

G 

840 

3.2 

10 

No RR 

— 

— 

610 

2.8 

10 

GW1 

SQ 

W 

6060 

0.0 

10 

GW1 

SQ 

E 

760 

1.9 

10 

GW1 

SQ 

G 

690 

2.4 

10 

GW1 

M 

W 

800 

0.0 

10 

GW1 

M 

E 

730 

2.0 

10 

GW1 

M 

G 

730 

2.5 

10 

GW2 

M 

G 

730 

1.2 

10 

GW3 

M 

G 

740 

1.1 

1 Stellar mass in Mq . 


2 GW approximations: GW 1 

Gair et al. 

2006 

, GW2 

Peters 

1964 

GW3 

Hopman & Alexander 

(2006a 







3 Coherence time: M = Mass prec., SQ = Self-quenching. 

4 Noise model: W = White, E = Exponential, G = Gaussian. 

5 Event rates in units of 10 - 6 yr — 1 . 


by the parameters of the host galaxy, in particular the MBH 
mass. Here we adopt a simplified one-parameter sequence 
of galactic nuclei, where the free parameter is M,, which to¬ 
gether with several additional fixed parameters define the se¬ 
quence. The M # -scaling is based on the empirical A/./rr rela¬ 
tion M. = Mo(a /(To)' 3 where a is the stellar velocity disper¬ 
sion outside the MBH radius of influence ry, = rj^GM^/a 2 , 
which encloses a stellar mass of order the mass of the MBH 
N+(rh) = ^hQ- The power law parameter /3 ~ 4 — 5 and the 


normalization Mg/a 1 ^ are determined empirically (Ferrarese 
|& Merritt] [2000 Ge bhardt et al.||2000l). It then fo lows that 


r h = vn (M./M o y~' A/p GM 0 /at$ 

Using this parameterization, and tb 


Alexander 


2011 ). 


tie approximation that the 
steady-state distribution is given by a BW76 cusp, the to¬ 
tal plunge and inspiral rates can be estimated from Eqs. 
and ( [32] ), 


FT 


10 


^3/2 2-kGMq \Mq 


M. 

M 0 
log Q 


3/0-1 


and 


JD t' 
-Lb: 


-Agw 


log 


io ,T 


^(M 0 /M.) 1/f} (c/4cr 0 ) 


(43) 


, ? 3/2 2ttGM 0 


(' Qo ) 


3/0—1 


(logQ) 1/s 


log VVhAcw (m/i log Q) ~ 2 / 5 Q 0 1/l3 °/ ( 4cr o) 


(44) 


where Q 0 = M,/M 0 , a G w = A G w (Vh log Q) -4/5 and 

A G w is a numerical factor which depend on the GW dissipa¬ 
tion approximation (Appendix [A). 

In our MC simulations, we adopted for simplicity /? = 4, 
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Fig. 17.— The total Plunge and inspiral rates as a function of MBH mass. 
The MC simulations (circles) agree wi th the ana lytic approximations for dy¬ 
namics without RR (solid lines), Eqs. <(46j, EZJ- Simulations with RR show 
that the contribution of RR is small: the discrepancy between the rates with 
and without RR does not exceed ~ 30% over 5 orders of magnitude in M # . 

/x/j = 2,rjh = 1. Mq = 5.4 x 10 6 Mq and cr o = 100 km s -1 . 

Thus 

r h =2(M. tMW ) 1/2 pc, (45) 

where M,^mw = M ,/4 x 1O 6 M 0 is the MBH mass scaled 
to the mass of the Galactic MBH. The rates as function of the 
MBH mass M. and the mass ratio Q = AI, /A/, are then 

<‘ = 3x 10- 4 M"^ 

_ logQ _ -i 

6.70- 0.25 log (M.,Mw) y ’ 

(46) 


and 

i^«5xlO- 6 M-^ 

_ (logQ) 1/5 _ 

4.65 - 0.25 log ( M.,mw ) - 2 log (log Q) /5 ’ 

(47) 

where we used the value Aqw ~ 0.029, corresponding to 
the GW dissipation approximation of jGair et al.|(|2006jl (Ap¬ 
pendix [A]). As shown in Figure [17] these analytic approxima¬ 
tions are in agreement with the results of the MC simulations 
over several orders of magnitude of M,. 

6. DISCUSSION AND SUMMARY 

The determination of the steady-state of galactic nuclei is 
a fundamental open question in stellar dynamics, with many 
implications and ramifications, and has been the focus of nu¬ 
merous numerical and analytical studies. In particular, current 
estimates of loss-rates vary over several orders of magnitude 
due to theoretical and empirical uncertainties. Previous stud¬ 
ies either used post-Newtonian A-body simulations, which 
are limited to small-A, or did not include the relevant rela¬ 
tivistic physics (Section [TJ. Building on recent progress in 
the formal description of RR as a correlated diffusion process 
(the 77 -formalism, |Bar-Or & Alexander|2014) l, we obtain here 
a MC procedure and analytic expressions for the steady-state 
distribution and loss-rates in galactic nuclei, taking into ac¬ 
count two-body relaxation, RR, mass precession and the GR 


effects of in-plane precession and GW emission. By cross¬ 
validating the analytic estimates and the MC results with a 
high degree of accuracy, and without the introduction of any 
free fit parameters, we are able to confirm our analysis and 
interpretation of the dynamics of the loss-cone in the context 
of our underlying assumptions. 

6.1. Discussion of main results 
The advantage of modeling RR by the 77 -formalism, over 


previous attempts by other approaches (Rauch & Tremaine 

1996; Hopman & Alexander 2006a; Giirkan & Hopman 2007| 

Madigan et al. |2011Merritt et al.|201I 

Antonini & Merritt 

2013'; Hamers et al. 2014; Merritt 2015a 

3 ), is that it allows to 


order relativistic 3D Hamiltonian. The resulting effective 
DCs, which are thus derived from first principles, are then 
guarantied to obey the fundamental fluctuation-dissipation re¬ 
lation and the cor rect 3D maximal entropy solution ( |Binney| 
|& Tremaine|2008| Section 7.4.3; Appendix|E]). 

These constraints on the functional form of valid DCs are 
critical, since the correct steady-state is the result of a fragile 
near-cancellation of two large opposing currents (the diffu¬ 
sion and drift); even small deviations from this relation (e.g., 
due to approximations, empirical fits, or reduction to lower 
dimensions), will result in large errors. For example, Hamers 
[et al.| ( [20T4| obtained the RR DCs from numerical simulations 
using an assumed functional form, D 1:j oc y 1 — j 2 , based on 
the fit of Giirkan & Hopman (2007 P and on the ad hoc expres¬ 
sion Dj j~ L Vjj, which is inconsistent with the fluctuation- 
dissipation relation and therefore leads to invalid steady-sate 
solution. This was then partially remedied by Merritt|(|2015ai 
who treated separately the Newtonian ( j — > 1) and relativis¬ 
tic^ —> 0) regimes. In the Newtonian regime, the [Hamers] 
et al.] p014| data was re-fitted to DCs that effectively satisfy 
the fluctuation-dissipation relation, which means that in the 
absence of a loss-cone, the dynamics asymptote to the max¬ 
imal entropy limit n(j ) = 2 j. However, in the relativistic 
limit j —> 1 , where the simulation statistics are poorer due 
to the smaller phase-space volume, Merritt] ( |2Q15aj) u sed an¬ 
alytic DCs based on the Hamiltonian model of Merrit t et al.| 
( 20111 ), which represented the stochastic background by an ad 
hoc dipole pseudo-potential and a recipe for switching its di¬ 
rection every coherence time. This recipe corresponded to the 
77 -formalism’s “Steps” or “Exponential ACF” noise (depend¬ 
ing on the exact switching procedure), which both converge 
to the same form in the j — > 0 limit (Bar-Or & Alexander! 
|2014[ Figure 1). As shown b y|Bar-Or & Alexander| ( |2014| Eq7 
42), in that limit I) n ss j 4 /T c and Dj ss (5/2 )j 3 /T c , where 
T c = 0 .5T c Uq R ( j = 1) /v'j ( j = 0). This indeed satisfies the 
fluctuation-dissipation relation, as any Hamiltonian model is 
guaranteed to do. These DCs are different from the ones de¬ 
rived by Merritt (2015a 1 , Djj oc j A /T c and Dj = 2 Djj/j, 


who implicitly forced the solution to 2D in-plane motion by 


setting sin i = 1 in the derivation (Merritt| [2015a Eqs.C. 8 - 
C.9). Therefore, these derived DCs satisfy the 2ITfluctuation- 
dissipation relation 2 Dj = dDjj/dj, rather than the correct 
3D one, 2 j D ;) = dj I), ,/dj. These DCs therefore imply the 
steady-state solution n(j) = const in the relativistic regime 
(assuming no loss-cone); this is not the correct solution for 3D 

9 We obtain a more accurate expression for Djj ( Appcndix[t)|, which fits 
torques measured in static wires simulations very well, over the entire range 
16 [ 0 , 1 ]. 
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orbital motion (the correct one is n(j) = 2 j). We note that 
the concatenation of two diffusion solutions, a 3D one for the 
Newtonian regime and a 2D one for the relativistic regime, 
may create an artificial discontinuity in the dynamical behav¬ 
ior at the interface. 

We have shown that the representation of stochastic dynam¬ 
ics near a MBH in terms of the streamlines of the probability 
flow provides a powerful tool for analyzing the loss fluxes, 
and leads to the identification of the exact separatrix between 
plunges and inspirals. We show that the typical sma of this 
separatrix, acw , yields an excellent analytical estimate for 
the inspiral rates found in our MC simulations (Figure [T7] >. 
This remedies the ambiguity in the identification of the crit¬ 
ical sma, which was used to estimate rates in previous stud¬ 
ies. We also explored the effect of different GW dissipation 
approximations on acw and the resulting GW inspiral rate, 
and found that the rates are robust to within a factor 2. Nev¬ 
ertheless, it is worth noting that the more accurate method 
of Gair et al. (|2006jl predicts EMRI rates that are more than 
twice higher than those of the commonly used approximation 
of |Peters & Mathews | ( | 1963] ). 

We have shown that GR precession plays a critical role in 
the dynamics of the loss-cone by efficiently quenching the RR 
torques. Conversely, in its absence all stars would rapidly 
plunge into the MBH, creating a depleted central cavity (cf 
Figure [l6|. This implies that GR precession is important even 
in systems that are effectively Newtonian, where at any given 
time only a very small fraction of the stars are on relativistic 
orbits. In particular, JV-body simulations of stellar dynamics 
and stellar populations near MBHs should include GR preces¬ 
sion even if the questions of interest are in the non-relativistic 
regime, so that plunges do not compete, or limit the lifetime 
of the stars. It is worth noting that GR precession has not 
yet been tested empirically for relativistic parameters larger 
than T = 2GM. /err > 8 x 10 6 (the double pulsar system 
PSR J0737-3039, Lyne et al.|2004j), whereas the S-stars in the 
Galactic center, for example, are already observe d to reach 
T = 2 GM,/c 2 r > 1 x 10 -3 atperiapse (star S14, Gillessen 


|et al.||2009a[ see also |Alexander||2006)l. Moreover, near the 
relativistic loss-cone (r p = Ar g , IGair et al. 112006b, where 
T —> 1/2, there is to d ate no empirical confirmation of GR 
precession (fWilIj|2006jl. The existence, dynamics and loss- 
rates of stars on such relativistic orbits can therefore probe 
GR precession in the strong field limit. 

We have shown that the influence of RR on steady-state 
loss-cone dynamics of compact objects is a small 0(1) effect, 
since the loss-lines for both direct plunge and GW inspiral lie 
well outside the region where RR is effective (e.g., Figuref!}. 
RR does introduce a correction to the steady-state distribu¬ 
tion and the loss-rates, which is at present small in compari¬ 
son to the astrophysical uncertainties, such as stellar density 
and mass function, the M,/a relation, and deviations from 
spherical symmetry. Eqs. ( [43] - [47] ) provide useful analytic es¬ 
timates for the plunge and tnspiral rates per galaxy, based on 
the NR-only approximation. The RR correction to the rates 
can be obtained by numerical integration of Eq. m Using 
a suite of MC simulations, we have verified in the context of 
our assumptions that the rate estimates are robust under differ¬ 
ent assumptions about the properties of the stellar background 
noise (smoothness, coherence time). This is in large measure 
a reflection of the limited role of RR in the presence of con¬ 
tinuous noise (jBar-Or & Alexander|j2014| Figure 1), which 
restricts the domain where RR is effective, so that slow NR 



log 10 (j) 



log a 


Fig. 18.— The suppression of the phase-space density of icy planetesimals 
around a MBH due to RR-driven interactions with a circumnuclear accretion 
disk. Top: The phase-space density and the disk interaction loss-line. The 
black contours delineate the region where RR dominates f D/f ! *- / = 

1.10). Bottom: The resulting steady-state a-distribution of the planetesi¬ 
mals for the NR-only case, and for the full dynamical model (top panel). 
The shaded regions denote the a-intervals where the loss-line crosses regions 
where RR dominates. 


remains the bottleneck and sets the rates. 

RR can substantially affect processes whose loss-lines in¬ 
tersect the phase-space region where RR dominates (D^ 
D RR ). The loss lines for tidal disruption of extended objects 
by the MBH, such as red giants or binaries, do lie closer to the 
RR line. However, it can be readily shown that neither class of 
objects is long-lived enough for RR-driven tidal destruction to 
play a dominant role in Milky Way-like galactic centers. Red 
giants are relatively short-lived, and the more extended and 
tidally-susceptible they are, the shorter their lifespan. Soft 
stellar binaries are destroyed by 3-body ionization before they 
are affected by RR-driven tidal separation ([Alexander & Pfuhll 

l2oni . 

One class of processes where RR may have more than an 
0 (1) effect is the hydrodynamical destruction (or removal) 
of objects by interaction with a large circumnuclear accretion 
disk. To demonstrate this point, we consider here as an ide¬ 
alized simple example the case of a massive accretion disk 
of radius = 2000r s (Goodman 20031 and a population 
of long-lived icy planetesimaTs around it, which are destroyed 
by several consecutive disk crossings (it is assumed that the 
number of crossings for destruction is N cross P <C Tnr)- In 
that case, the critical angular momentum for destruction is 
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jd, = \J2GM % Rd = 15.8 ji c , which is large enough to inter¬ 
sect the RR-dominated zone (Figure [18] top). As a result, the 
differential sma distribution of the planetesimals is depleted 
below a < 100 mpc and strongly so below a < 10 mpc (Fig¬ 
ure [18]bottom). 

Although RR is typically inefficient in driving stars all 
the way to the loss-cone, it can randomize orbits in the 
phase-space regions where it dominates, even when the NR 
timescale is longer than the system’s age or the lifespan of 
the stars. This may be a key element in solving the “paradox 
of youth” (|Ghez et al.||2003Jl in the Galactic Center, where 
young B stars are observed on tight orbits around SgrA* (the 
so-called “S-star” cluster). The leading formation or migra¬ 
tion scenarios for the S-stars predict non-isotropic initial ec¬ 
centricities (see reviews by |Alexander|2005] |2011) >. However, 
most of the S-stars are in the RR-dominated phase-space re¬ 
gion (Figure [TJ, and so substantial evolution and isotropiza- 
tion of the initial eccentricities is possible. However, many 
of the S-stars are short-lived, and some are close in phase- 
space to the Al-suppressed region. A detailed analysis using 
the ^-formalism, which can treat this intermediate dynamical 
regime rigorously, has yet to be carried out. 

6 . 2 . Limitations and caveats 

The applicability and validity of our results are limited by 
several simplifying assumptions. We assume a non-spinning 
MBH surrounded by an isotropic (on average) non-rotating 
single Keplerian power-law cusp of single-mass stars. We as¬ 
sume the dynamics are dominated by single star interactions. 


that is, we neglect the possible effects of binaries, or the con¬ 
tribution of non-stellar massive perturbers , suc h as gas clouds, 
cluster or intermediate mass BHs (Perets et a l.[2007} . 

In our MC simulations we used NR and RR DCs, which are 
based on a fixed isotropic, single power-law BW76 model. 
These DCs are therefore not self-consistent with the steady- 
state solution. However, we showed that the BW76 cusp is 
a good approximation (within a factor of two) for the steady- 
state solution in the relevan t reg ion (down to asw < clgw for 
M. < 10 9 M.; see SectionpAk. In addition, it can be shown 
that for the RR DCs the isotropic fluctuation-dissipation rela¬ 
tion holds even for the non-isotropic case as long as the total 


angular momentum of the system is zero (Bar-Or & Alexan 
|der||2014J). Since the fluctuation-dissipation relation results" 
from the symmetries of the Hamiltonian, it is reasonable to as¬ 
sume that the same will hold also for the NR DCs. This means 
that small non-isotropies will only result in small magnitude 
changes of the flux and have little effect on the steady-state 
distribution. 

Finally, the RR DCs are based on simplified (single 
timescale) background noise models that can be treated an¬ 
alytically, and simple coherence models that are functions of 
the sma only. However, we were able to show using MC sim¬ 
ulations that the results are largely independent of the exact 
noise model as long as it is continuous (i.e., not white noise). 


T.A. acknowledges support by the I-CORE Program of the 
PBC and ISF (Center No. 1829/12). 


APPENDIX 

A. THE ANALYTIC GRAVITATIONAL WAVE LINE 


The phase-space of the relativistic loss-cone is divided by a separatrix into an outer region where streamlines end as plunges, 
and an inner region where strea mline s end as inspirals (Figure [4]. Since the probability current along an infinitesimal bundle of 
streamlines is conserved (Section|3.5}, it can be evaluated at any convenient point along the flow, in particular at J > J) r , where 
GWs are negligible, and the streamlines are identical to those in the absence of GWs. Therefore, the inspiral rate is estimated by 
locating the no-GW plunge streamline that corresponds to the separatrix, and identifying the sma of its terminal point (a p , ji c ) as 
the critical (maximal) sma for inspiral, acw- The inspiral rate is then simply the integrated plunge rate in the absence of GW, up 
to a Gw . i.e., Rl ot = R™ GW {a GW ). 

The separatrix streamline, dE/dJ = Se/S j, can be evaluated by noting that the flow in the ^-direction is mostly due to GW 
dissipation, 

S E = n (E, J) E/taw , (Al) 


where tow = a/ da tv is the GW orbital decay timescale. In the presence of GW, the innermost plunge streamline initially 
approaches ji c from j ~ 1 at nearly constant a, and then turns over to lower a and runs nearly parallel to the loss-line ji c (a) = 
4^/rg/a, until it terminates when j -a ji c at some small enough terminal sma, ai nsp , where the orbital motion is effectively 
circular (i.e., inspiral, cf Figure |4jl. For example, for parabolic orbits, the minimal possible value for a lnsp is obtained where 
ji c = 1, i.e., for ainsp = 16 r g . As we demonstrate below (Eq. <|A9}), the separatrix solution is a function of Gq nsp only via 
ctinsp/ftmax <C 1 , and is therefore independent of the exact choice of ai nsp . 

The flow in the J-direction is approximately constant in J with a typical timescale Tj (Eq. (| 1 6}), 


n(E) 

Tj log {J/Jic) ' 


(A2) 


Since RR is negligible near j = ji c , Tj can be approximated as Tj (a) ss 9 1 Q 2 P (a) / (N (a) log Q) (see Eq. ([19}), where 9 1 is 
a numeric pre-factor that depends on the cusp density profile (0~, ss 1/8 for a 7 = 7/4 (BW76) cusp). 

It then follows that separatrix is the solution of the differential equation 


dE 

dJ 


Se 




E 

tow 


log 



(A3) 


with the boundary condition J(Ei nsp ) = Ji c (where sp = GM,/2ai nsp ). 

An exact expression for the GW dissipation can be obtained only numerically. Some useful analytical approximations are 
available (see Gair et al. (|2006|) for a comparison between the different techniques). The simplest expression was obtained 
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by Peters & Mathews | ((1963), who assumed a point-mass objects moving on a Keplerian orbit. In this approximation the GW 
timescale is 


taw = tt - 


1 5 Q r„ 


2r, 


7/2 


P(a),f(e) = 

T, 


\a G w\ 27t 64 /(e) a 

r c , the periapse r p is rela 
(Eq. (|A3)) can be written in terms of x = a/a max and s = J/Ji c , 


l+e 


-7/2 


73 


37 


1+ 24 6 + 9lT 

r 2 / 72 


(A4) 


In the limit J -A J c , the periapse r p is related to the angular momentum by r p /r g = 8J 2 /J 2 C and the streamline equation 


= AdX 1 2 s 6 log (s) , 
as 


where the competition between the GW dissipation and NR diffusion is expressed by the parameter 

The boundary condition at the terminal sma is 


X\ — X (s — 1) — dinsp /ttmax ^ I 5 

and the plunging branch of the solution is give by 


•(s) = 


3-7 , A d (3 — 7 ) (s 5 — 5 log s — 1 ) 
'• +-25^- 


1 1/(3 — 7 ) 


(A5) 


(A 6 ) 


(A7) 


(A 8 ) 


In the phase-space region where GW is negligible, the streamline defined by Eq. ( |A3) is approximately constant in E. There¬ 
fore, dew can be identified by the sma at s /g> 1 (j —>• 1 ), that is 


OGW ~ ®max 


lim x = (x\ 7 + ^ —Ad\ 

S— >00 \ Zb J 


l/(3-7) 


For a steady-state cusp (i.e., BW76, 7 = 7/4), x\ 7 <C (3 — 7 )^ 73 / 25 , and therefore 

-4/5 

, / 2 V h iUft ^ 1 

O.GW 




(A9) 


(A10) 


where the value of Aq W depends on the specific GW dissipation approximation that is assumed. For the Keplerian approxima¬ 
tion A% w = (7t/(1)/3200) 4 / 5 ~ 0.013. 

As shown by Gair et al. (2006), the Peters & Mathew s](|1963[ ) estimate can be improved by using a “semi-relativistic” approxi¬ 
mation, that is using the fully r elativistic orbit in place of the K eplerian one in Peters & Mathews (|1963) equations. In the limit 
e —s- 1, this ap proach, used by | Hopman & Alexander (2005 1 , amounts to re placing the Keplerian r p /r g in Eq. ( |A4) with the 
relativistic one (|Gair et al.|2006[~ 


r p/ r g = 4S + \/s 2 - 1 'j . 

Thus the GW separatrix (Eq. ( |A5) ) is replaced by its relativistic version, 

dx .7/0 . o / . r -7 ; x\ 7/2 

ds 


2 ,//2 A c X 7 2 ^1 + y/l — 1/s 2 ^ S 6 log(s) 


(All) 


(A12) 


where a G w can be solved numerically. As shown in Figure |T9l agw can be approximated for 7 = 7/4 by Eq. ( |A10| > with 
Aq W = 0.022 and for the more accurate treatment of Gair et al. (2006), A G w = 0.029, which is adopted in this study. 

B. MONTE-CARLO SIMULATIONS OF LOSS-CONE DYNAMICS 

We summarize here the assumptions underlying our Monte-Carlo (MC) simulations, the details of the implementation and the 
derivation of the steady-state configuration and loss-rates. 


B.l. Assumptions and procedure 

The MC simulations in (a,j) assume a fixed background model, whose properties define the NR diffusion coefficients, an 
RR coherence timescale model and a background noise model. The stellar background is approximated as a power-law cusp 
with enclosed number of stars N(a) = -A max (a/a max ) 3 ~“, extending between a m i n and a max . The accessible phase-space for 
the test stars extends between a reflecting boundary at a out = a lnax (evaporation), and an absorbing boundary at «in < amin 
(destruction). The reflective boundary at a ou t ensures that the long-term distribution of test stars converges to an isotropic 
distribution in j, n(j) -A 2 j, which is the assumed distribution of stars far from the MBH (note that this is not guaranteed for an 
absorbing outer boundary, even when the test stars are introduced into the simulation isotropically). The extension of phase-space 
to small oin allows inspiral trajectories to be tracked down to a tight enough sma where their ultimate fate as EMRIs is certain 
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Fig. 19.— The semi-major axis o,gw °f the outermost inspiraling streamline for a BW 76 cusp (7 = 7/4) for different MBH to star mass ratios and for 
different approximations of the GW dis sipa tion. The approximate expression (Eq. jAlOl ) for ac; 14 / (lines) is fitted to values of (ic; w found by numerically 
integrating the streamline equation (Eq. lA3j) (circles). 


(cf Figure [T5j ) phase-space extends in j between a reflecting boundary at j = 1 and an absorbing boundary at the last stable 
orbit (LSO), jr.so* defined by the locus Jlso{o) = 7 / a/r g jLSO = 4r g c (the LSO for a zero energy orbit in the Newtonian 
approximation; in the limit e —>• 1 it corresponds to a critical periapse r^so = «(1 — e) = 8 r g ). 

An MC run starts by injecting a test star in some initial phase-space position (ao, jo). This is typically chosen randomly, 
either just below a max , with j distributed isotropically above j [.so* to simulate a star diffusing from an isotropic galaxy into the 
MBH’s radius of influence, or isotropically in the cusp’s bulk, according to the power-law cusp distribution, to simulate the initial 
conditions of an iV-body simulation (see below), or the distribution of tracer stars (e.g., red giants), which mirror the distribution 
of cusp stars. The star is then evolved in small time increments d t, taking into account the stochastic changes in energy and 
angular momentum due to NR, the deterministic changes due to GW dissipation, and the random changes in angular momentum 
due to RR. The star is tracked in phase-space, and its position is recorded by snapshots taken at fixed intervals At dt. 
Ultimately, after surviving for time t s , the star leaves the system at some terminal phase-space position (asi, j 1 ) as a result of one 
of four possible outcomes. (1) Evaporation. The star’s sma crosses a ou t to a larger sma. This happens for the majority of test 


stars. (2) Inspiral. The star’s orbit decays until it crosses a». 
directly, at some a p > a, insp . 


= eomin mpc (e <C 1). (3) Plunge. The star crosses the LSO 
(4) Finite lifespan exceeded. This is r elev ant for stars, which are limited by stellar evolution or for 
binaries, which are also limited by dynamical evaporation (Section l 6 Tj), but not for compact remnants. Note that the branching 
ratio between plunges and inspirals is independent of the exact value of the sma chosen to distinguish between the two outcomes 
(parameterized by e), since their respective terminal phase-space positions are clearly separated (min a p i ung e A> ftinspirai, cf 
Figure [j~5) e = 0.1 was typically used here). 

Once a star exits the system, a new star is injected (reflection at a ou t is equivalent to injectio n of a new star at (a ou t, j \ )). 
This is repeated over some long accumulated time T s ; m 2 > t s (typically T s ; m ~ 1007).;, Section 13.41). All test stars simulated 
over T s i m , apart for the last one, which is omitted from the analysis, reach a definite outcome. For each test star we record its 
trajectory in phase-space (in coarse At resolution); the total time it spent evolving in phase-space, t s \ the nature of the final 
outcome (evaporation, inspiral, plunge or end of lifetime); and its initial and terminal phase-space positions. The procedure is 
repeated as needed (typically 10 3 — 10 4 x T s i m ), until enough test star statistics are collected. The snapshots of the phase-space 
positions, the survival times and the final outcomes are then used to estimate the steady-state configuration and loss-rates, as 
detailed below. 


B.2. Representation of physical processes 
The form of the RR DCs depend on the background noise model and the precession of the test star (|Bar-Or & Alexander! 


2014). We assume three optional noise models: white noise (equivalent to no precession), C° noise with exponential ACF (an 
Ornstein-Uhlenbeck process), and smooth C°° noise with a Gaussian ACF. Prograde GR in-plane precession is modeled by the 
1PN approximation isgr = 3 v r {a){r g /a)/j 2 , where v r (a) is the Keplerian radial (orbital) frequency. Mass precession is given 


exactly for an a = 2 cusp, j) = — \N(a)/Q\v r j/(1 + j) (Merritt et ak 20111, and for other values of a by polynomial 

approximations of the exact integral (Alexander |2005| ). The magnitude of the RR DCs reflects the strength of the RR torques, 
t n ~ 0.28 y/T^jy/NpajGM^/a, as derived from static wire simulations (Appendix IdI). 

The RR DCs have explicit analytic forms that are easy to evaluate. However, the NR DCs (Appendix O involve multiple 
integrations that are too computationally expensive to perform on the fly. We therefore calculate them exactly beforehand on a 
25 x 25 evenly-spaced logarithmic grid extending between a m i n and a max and ;j ui \ u = 10 ~ 3 and j max = 1, and then bi-linearly 
interpolated to any (a. j) as needed. 

Three optional analytic perturbative es tima tes for the rate of GW dissipation of energy and angular momentum were studied: 
those of Peters] ( j 1964| ), |Hopman & Alexander] ( j2006aj ), and Gair et ak ( j2006| . 
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Each MC step, the time-step was chosen as A t = min (AtNR , Ataxy), where the time-steps for the NR and RR were 

chosen by criteria similar to those used by Shapiro & Marchant] ( jl 978| for NR, and the time-step for GW was chosen to be a 
small fraction of the GW dissipation times mm(E/E qw , J/Jgw)- 
Two optional approximate background coherence timescale models are considered, which are assumed to be functions of a 
only: /-averaged mass precession, Tm = \/ k/ 2v^} (2a, \/l/2), and self-quenching Tsq = Qv~ 1 {a) /A(2a), where the 
approximate numeric pre-factors (Bar-Or & Alexander 2014, footnote 7) were evaluated here specifically for a = 7/4, but are 
generally insensitive to the exact value of a. Note that GR precession is not included in background coherence models of the 
MC, but it is included in approximate form (;/ prec = vm + vgr\) in the analytic modeling of the coherence time (Section^. 


B.3. Steady state rate estimates 

The interpretation of the MC results in terms of loss-rates depends on whether the test stars represent the underlying background 
cluster that is generating the NR and RR perturbations, or whether they are a separate trace population that is affected by the 
perturbations, but does not contribute to them. 


B. 3.1. Test stars as background 

In statistical steady-state, stars that exit the system (a;,, < a < a ou t) are replaced at a rate that keeps their time-averaged 
number N fixed. A star that evaporates from the system back to the infinite reservoir at a > a ou t (reflection at a ou t is treated 
as evaporation), is replaced by another star from the reservoir, so there is no net current through a out due to evaporation. The 
situation is different for stars that end up in the MBH, whether by plunge or inspiral. Since they are permanently removed from 
the system, a net current of stars through from the reservoir into the system, is required to compensate for their loss. Both 
the stars that evaporate and those that are losj^j contribute to the total mean number of stars in steady-state. In our models, this 
number is an input parameter, determined by tEe assumed background cusp. 

Designate by P& the probability (branching ratio) for outcome k = 0,1, 2,..., where k = 0 denotes evaporation and k > 0 
denote the various loss channels, so that N /i; P k = 1. The MC statistics provide estimates of the branching ratios, = n;./n s ; m , 
where ?r s ; m is the total number of test stars whose phase-space trajectories were simulated, and ri^ the number of times outcome k 
has been reached. This translates to event rates by requiring that the total number of stars be on average A = A r / :: = T ^tk, 

where T/c is the rate of outcome k, and tk = t\ l JA is the mean survival time in the rik simulations that had outcome k. 
The rate for each channel is related to the total rate of all outcomes, T, by T^ = TP/,., where T = N/t s and t s is the overall mean 
survival time, irrespective of outcom It then follows that the event rates are 

r fc = (N/i s )Pk . (Bi) 

The total replenishment rate that is required to keep the system in steady-state is then the sum over all the loss channels, 

Tloss =J2 Tk - (» 2 ) 

fc >0 

Loss-rates estimated by direct A'-body simulations (Merritt et al.] |2011) [Brem et al.| |2014| ) can in principle be compared to 
the rates derived by MC simulations of scaled-down nuclei, where the test stars are treated as representative of the background. 
However, this comparison is complicated by the fact that the A-body systems are not necessarily in steady-state (the initial 
configuration may not be the steady-state one, and / or stars lost in the course of the simu lation are not replenished), and do not 
have fixed boundaries or boundary conditions. The comparisons discussed in Section |4~2| and Tableware approximate. The MC 
loss-rates were estimated for a non-equilibrium cusp that corresponds to the initial conditions of the A-body simulations, and the 
MC rates were compared to the rates early in the simulations, at times where the A-body configuration is still close to its initial 
state. 


BAA. Test stars as a trace population 

It is of interest to consider how a small population of tracer stars, which are injected into the cusp by some dynamical or evolu¬ 
tionary mechanism, evolves dynamically with time, and is lost via the various channels. Some possible injection mechanisms are 
capture by tidal separation of an incoming binary (jHills[l988j), in which case the injection point in phase-space is a tight eccentric 
orbit deep inside the cusp; the formation of a red giant when a background star evolves off the main sequence, in which case the 
injection point reflects the background distribution of progenitors; or the formation of massive blue giant in a fragmenting gas 
disk, in which case the injection point is a low-eccentricity orbit. 

When the test stars in the MC simulation represent a tracer population, the total rate of all outcomes is determined by the 
assumed injection rate, and is no longer related to the total number of background stars. In this case the MC does not predict the 
event rates Tfc, but rather the branching ratios Pk. 


10 To simplify bookkeeping, a test star that would have wandered back and 
forth across a ou t and is finally destroyed by the MBH is not counted as a 
single star (i.e., a single trajectory). Once it leaves the cusp (a > a ou t), it 
is considered as evaporated. The next crossing into the cusp a < a 0 ut is 
identified as the beginning of the phase-space trajectory of a new star. 


11 This follows from summation over all channels: N = A,. I 'k.th 
T A h. Hfcffc = r Afc( n fc/ n sim)( n fc Ay ^7// = r(Afc,y )/ n sim 

rt s . 
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C. TWO-BODY DIFFUSION COEFFICIENTS 

We derive here the relations between DCs in velocity space and the DCs in angular momentum and energy space for stars 
orbiting in a spherical potential. 

Consider a star of mass m moving in a spherical potential, 6 = 6 (r), with velocity v. The binding energy and angular 


momentum are 

E = <f)(r) — ^v 2 , (Cl) 

J = r x v . (C2) 

Note that here E is the positively defined orbital energy and </> (r) is positively defined potential. Due to gravitational encounters 
with the field stars, the star changes its velocity, v, to v' = v + Av. Consider the orthonormal basis 

v = v /v , (C3) 

J = J/J = rxv/|rxv| , (C4) 

w = vxJ/|vxJ| = v x (r x v) = (nr — v r v) jvt . (C5) 


In this basis, the change in velocity is 
where 

and 

The change in energy is 

A £=-i(>/ 2 -» 2 ) 


Av = Aiimv + Avj_, 
Avj_ = AnjJ -f Av w w , 

Anj_ = yj (A Vjj 2 + ( Av w ) 2 . 


( Av f - v • Av = (Any) 2 - ^ (An) 2 - nAnj . 


(C6) 

(C7) 

(C8) 


(C9) 


The position vector, r, is 

r = r ( v t /v) w + r (v r /v) v , 

where v r and v t are the radial and transversal velocities. Therefore, the change in the radial velocity is 

. Av • r v r . v t A 

Av r = -= —A tin H- Av w . 

r v v 

The change in the transverse velocity up to second order in Av jv is 


Aug Av u 

A Vt = Vt - V r - 

V V 


1 A v 2 j 

-\ - - 

2 v t 


The change in the angular momentum is 

A J = r x Av = J 


Ann 


—An,,, j J T J 


Avj ( v 


v t 


and the change in angular momentum magnitude (up to second order in Av/v ) is 


V V 


1 r 2 2 
--An 2 ,. 

2 J J 


Using Eqs. (|C9[> and (|C 1 4[>. we can obtain the local (orbital phase dependant) DCs in terms of the velocity DCs, 


(A E) = -^ <((An||) 2 ^) - i ^(Anj_) 2 ^ — v (Any) , 
((Af?) 2 }=n 2 ((Ann) 2 ), 

< AJ > = ^ < AW H) + ij(( Aw± ) 2 ) ’ 

(( A ^) 2 ) = ^ <( A -ll) 2 ) + \ ( r2 (^) 2 ) 

(Ai?AJ) = -j((An||) 2 ) , 


(CIO) 

(C11) 

(C12) 

(C13) 

(C14) 

(C15) 

(C16) 

(C17) 

(C18) 

(C19) 
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; we omitted higher order terms in Av/v, and used (Av w ) = (A vj) = 0 and (Av 2 r ) = (Av 2 ,) = (Av 2 ) /2. 
' - - - -- — ' 8 ) 


The local velocity DCs are (Binney & Tremaine|2008] 


(A-l’ii ) = -K m + ma ( dv a ^f a (v a ) , 
'"'’a J 0 v 


TO, 


Avf\ ) = -K 




f v v 4 

/ dv a -jfa(v a ) + / dv a V a f a (v a ) 

Jo v Jv 

r /I ?; 2 ?; 4 \ Z 100 

J dv a ^ - -§J / a (fa) + 2 J dv a V a f a (v a ) 


(C20) 

(C21) 

(C22) 


where /c = (47rGm a ) In A, and where f a (v a ) and m a are the velocity DF and mass of the field stars, and where we assume that 
the velocity DF is isotropic. In that case the velocity DF can be written in terms of the orbital energies of the field stars. Using 

f a (v)vdv = -f(E)dE, (C23) 

we obtain 


(AE) = k 


(AE 2 ) = z - 


m f 

- / ( v a/ v ) fa (E a ) dE a f a (E a ) dE a 

m a J E J- oo 


KV 


/•(j) pE 

/ dE a ( Va /vf f a (E a ) + / dE a f a {E a ) 
J E J—oo 


(AJ) =K < (— + l') [ dE a (v a /v) fa(E a ) 

( v A \m a J J E 


and 


(. AEAJ) = -J-k 

o 


r dEa (v a /vf f a {E a )+ f dE a f a (E a ) 
J E J —oo 


-1 


The corresponding orbit-averaged DCs are given by 

D e = 2P 

D ee = 2 P 
Dj = 2 P- 1 
Dt.t = 2 P" 1 


D ej = 2 P 


-1 


(A E) dr/v r , 
((AJT) 2 ) dr/v r , 
(AJ) dr/v r , 
((AJ) 2 ) dr/ tv, 
(AEAJ) dr/v r . 


(C24) 


(C25) 


+ kf J dE 'f a ( Ea '> 3{v a /v) - (v a /vf + TfJ j dE a f a (Ea) 


(AJ 2 ) = 4 I J 2 [ [Va/vf dEafa (E a ) - J 2 ^ dE a (v a /v) f a (E a ) 

I J E J E 

+ 7 —^~ J dEa (s (v a /v) - (v a /v) 3 ^j f a (E a ) + ^r 2 V 2 J dEafa (£^a)| , 


(C26) 


(C27) 

(C28) 

(C29) 

(C30) 

(C31) 

(C32) 

(C33) 


where P is the orbital period. 
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Assuming that the potential is Keplerian, 0 = M,/r, the energy and angular momentum are 

M, 1 2 1 J 2 _ GM. 

- V “ 2 Vr ~ 2 Z 2a ’ 

J = rv t , 


(C34) 

(C35) 


where a is the sma. Using x = {r/a — 1) e 1 , the orbital average is 

1 r 1 i± w 

(D) 0 = - D^Z^dx, (C36) 

7T J_ x VI - X 2 

where e = yfl — J" 1 1 J 2 is the eccentricity of the orbits and J c = \/GM % a is the maximal (circular) angular momentum. The 
orbital-averaged DCs are therefore 


where 


De/E = —T 0 -|-Tno , 

m a 

, ~ 4 4^ 

Dee/E =-Ti 3 _i +-T 0 , 


Dj/J c = 


5-3 j 2 , 2 m + m a 1 

-To 1 0 — J —~- 1 111 + 1 310 — o 1 330 

12 2 m„ 3 


/j, 


Djj/J c — —j^-r 0 + -j 2 Ti3i — -j 2 r m + 2 r 3 i 0 — ^r 330 , 

Dej/ ( EJ c ) = —-j (To +Ti 30 ) , 


T 0 = k f dsfa(sE) 
J — OO 


(C37) 

(C38) 

(C39) 

(C40) 

(C41) 

(C42) 


and 


r ijk (e , j) = 2 


1+k—i ' 


r 2/(l+ex) 


dx 


dsfa ( sE) 


{r/a) 1 


1-1 


x/ZZ 2 {v 2 /E) f 


{v a /v) J 


(C43) 


We can write the Ty* functions in terms of the Cohn & Kulsrud (1978|> Fi functions (see the definitions after Eq. (24) there) 


F 0 = ET 0 , 

(C44) 

Fi =ETno , 

(C45) 

F '2 = ETni, 

(C46) 

F 3 = ET 3W , 

(C47) 

f a = et 13 _ 1 , 

(C48) 

E 5 = ET i 3 o , 

(C49) 

F 6 = ET i3i, 

(C50) 

Fr = ET 33 o . 

(C51) 


It is sometimes useful to consider the diffusion in the dimensionless normalized angular momentum, j. For a general coordinates 
transformation x' = x' (x) the new DCs are given by (e.g., Risken 1989) 


If = dX ’i D I 1 d2x 'i 
1 dxk k 2 dx r dxk 
, _ dx\ dx' m 


E r k 5 


dx r dxk 


(C52) 

(C53) 


Thus, in the E, j coordinates the ^-related DCs are 


Dj = Dj/J c + l -jD E /E + X -D E j/ {JcE) - \jD EE /E 2 , (C54) 

Djj = Djj/Jl + i j 2 D EE /E 2 + jD EJ / ( J C E) , (C55) 

Dej/E = Dej/ {J c E) + -jDee/E 2 . 


(C56) 
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Similarly, for the E, R 


j 2 coordinates the R -related DCs are (see also Cohn & Kulsrud 1978 Cohn|l979 


2 De , n D EJ 


Dr — “2-jDj + Djj — 2 jDj/J c + Djj/J c + j — —f 2 j 


E 


J,E ‘ 


Drr = 4 fDjj = Aj 2 Djj/J 2 c + f D EE IE 2 + 4j 3 
D E r = 2jD E j = 2j—^ + J 2 Dee/E 2 . 


(C57) 

(C58) 

(C59) 


D. THE STATISTICAL PROPERTIES OF THE RESONANT TORQUES 

We describe and measure here the residual torques acting on a test star due to the near-Keplerian orbits of the background stars. 
The torque vector r depends on the test star’s angular momentum J and argument of pericenter uj. We discuss the symmetries 
of these torques and the scaling with the test star’s orbital parameters, and empirically measure the torques by static random 
background simulations. 


D.l. Geometrical description 

Consider the angular momentum vector in some fixed reference frame 


J = J L, 


\J J 2 — J 2 COS (j) 

= ( V J2 - J l sin < 

J , 


where £ = (£ x ■ £ y ■ C) is the unit vector in the direction of J in a fixed Cartesian reference system (x. y. z). The torque is 
derived from the Hamiltonian by r = J = { J, H}, where {... } denotes the Poisson brackets. It is more convenient to represent 
the torques in an orthonormal spherical coordinate system (J, cf> , u), where u = l z = cos 6, with the associated unit vectors 
Bi = ( dJ/di) / \d3/di\ for i £ {J,(f>,u}. The change in the angular momentum’s magnitude is then tj = J, and the changes in 
its direction are described by the angular torques, r q j = J& and t„ = Ju. The transformation of the torque vector from spherical 
to Cartesian coordinates is given by 


t = Tjej 


+ T0\/l 


1t 2 <50 + 



(Dl) 


In the reference frame of the orbit we can define the orbital torques in the direction of the semi-major axis, r a , and the semi¬ 
minor axis, Tj,, 


r a = T ■ a = \/1 — it 2 sin LOTij, - t cos ujt u , 

Vi — u 2 

Tb = T ■ b = \Jl — u 2 cos ujTrf, H—sin ujt u , 

VI — u 2 

where a and b are the direction semi-major and semi-minor axes 

a = &<!, sin ui — e u cos ui , 
b = e<f, cos u> + e u sin w . 


(D2) 

(D3) 

(D4) 

(D5) 


D.2. Statistical properties 

Define the typical (Poisson) torque ( Giirkan & Hopman|2007j > 


TJV = 


y/Nj2a 


-GM* = J r v, 


Q 


The mean squared values of the torques are given by 

(tj)/?n=T\\ (a, J) , 

(-4) /t 2 n = ^ l T + ( a > J ) + T - («’ J ) cos 2w ] ) 

(r 2 ) /t% = ^ s/l — p 2 [T+ (a, J) - T_ (a, J) cos 2w] , 


(D6) 

(D7) 

(D8) 

(D9) 


and the cross terms are 


(tjtV = ( TjT u ) = 0 , 
(TuT^) /t% = T_ sin (2w) /2 


(DIO) 

(Dll) 
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Fig. 20.— The mean square of the torques (normalized to the typical (Poisson) torque f) as function of the angular momentum. The torques measured from 
static wires simulations (circles) are approximated by polynomial fits (solid lines). As expected from considerations of symmetry, for a circular orbit (J = J c ), 
there is no torque in the J direction and both perpendicular torques ( 77 ,, 77 ,) are equal, i.e., Ty = T_ = 0. As J —¥ 0, the orbit’s geometry approaches a rod, and 

the torque in the a direction vanishes because its lever arm goes to zero, while the torque in the J and b directions become equal, i.e., Tn = T_|_ = T_. 


where we defined 


T \\ = ( t j) Z^N ■■ 

(D12) 

T+ = (( T o) + ( T b )) AiV > 

(D13) 

T- = {( T b) ~ ( T a)) /^N ■ 

(D14) 


Note that since the torques (tj), (t^) and ) are measured in the orbital plane, they have no angular dependencies, and neither 
do Tj|, T + and TL. 

The symmetry of the orbits is such that for a circular orbit (i.e., J = J c ) there is no preferred direction in the orbital plane and 
therefore (r^) = (r^) and Ty (a, J c ) = T_ (a, J c ) = 0. For a radial orbits (J = 0), r a vanish and (r^) = therefore, 
T\\ (a, 0) = T + (a, 0) = T_ (a, 0). 

D.3. Measuring the torques 

We used static wires simulations (e.g., |Giirkan & Hopmanj|2007j) to measure the ^-dependence of the residual torque (j = 
VT^) on a test orbit with sma a and eccentricity e. This was carried out by simulating the background as many fixed 
Keplerian wire orbits, and measuring the three components of the orbital torques r a , Tb and tj, for many independent random 
realizations of the background, integrating over the orbit of the test star and the orbits of the field stars with the efficient iTouma] 
et al. (2009 1 algorithm. We decomposed the measured orbital torques to Tn, T + and T_ and fitted them to a second order 
polynomial in j. The best-fit results are (Figure [20]> 



T\\* 

0.08 (1 — j) (1 — j/8) , 

(D15) 



0.08 (1 — j) (1 + j/ 4) , 

(D16) 



0.08 (1 - 5j/8) . 

(D17) 

The residual torque in the J direction is therefore 




3 

/{t$)/t n W 0.28 VU - j) (1 - j/8) W 0.28 Vl - j ■ 

(D18) 

This analytic fit is consistent with the 

Giirkan & Hopman 

2007 1 results (fitted to yj(rj )/tn ~ 0.25\/l 

— j 2 see Eq. (13) there) 

but is better, since it matches their data over the entire 

j-range. For the out-of-plane torques we obtain 



Ao 

A 

III 

7 

r a + T b ~ 0.28\/l — 5j/8. 

(D19) 
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This is different from Giirkan & Hopman (2007) results (fitted to \J (r/)/Tv ~ 0.28(3/2 — j 2 ) see Eq. (14) there). However, 
the discrepancy can be traced to an error in their randomization procedure (Eqs. (9-11) there) which is not truly isotropic. 

Using A’-body simulations, Eilon et al. (2009) measured the isotropic averaged residual torques. They defined and measured 
the /-averaged quantities 


and 


a 2 = h ^ p2,Jl ) = 4,rl ■ 


(D20) 


(D21) 


Eilon et al. 


( |2009p 


Averaging the results here over j, we obtain /3 S ss 1.0 and j3 v « 1.7 for 7 = 7/4, in agreement with 

E. THE MAXIMAL ENTROPY PRINCIPLE 

The maximal entropy principle (MEP) has been shown to be a powerful tool in determining the steady-state (or quasi-steady- 
state) of dynamical systems. In particular, it was studied extensively in the context of collisionless self-gravitating systems (e.g., 
Lynden-Bell |1967| . Here we examine the more restricted problem of a near-Keplerian system. We prove that the RR DCs used 


in this study comply with the MEP. For NR, we prove that the J-only NR DCs comply with the MEP, since for stellar systems, 
the MEP is relevant only when the interactions conserve energy. 

In a system where a central object of mass M. dominates the potential (e.g., planetary systems, nuclear clusters), stars move on 
nearly Keplerian orbits. That is, the potential is almost regular and the orbital elements are almost constant. In particular, since 
the potential varies on much longer timescales than the orbital time, the potential can be orbit-averaged and the Keplerian energy 
is conserved. 

The entropy of the system is given by 


S = ~J d 3 rd 3 vf (r, v) log / (r, v) , 

where / is the stellar DF. The Keplerian energy distribution relative to the MBH, n(E), is conserved 


n(E;f) = J d 3 rd 3 vf (r ,v)6^E 


GM. 1 o 

— + v 

v 2 


(El) 


(E2) 


This implies the conservation of the total Keplerian energy E. [/] = f En (E; f) dE and of the total number of stars N tot [/] = 
f n(E- : /) dE. Conservation of the total energy E to t = E. [/] + 2?+ [/] then implies the conservation of the total potential energy 
due to the interactions with the stellar background. 


where 


E *{f} = \ J d 3 r J d 3 vf (r, v) ip (r; /), 

3 //(r'X) 


r — r' 


1/ (r; /) = GM* J d 3 v' J d 3 v 

is the star-star potential. The total angular momentum is also conserved, 

L fo t[/] = J d 3 rd 3 wf (r, v) L , 

where L = r x v. Using the Lagrange multipliers /3, b and A (E), we write the target function 

S = S + /3 (/«v/ (r' ; v') ip (r; /) - 2 E^ + b • (^J d 3 rd 3 vf (r, v) L - L tot 


(E3) 

(E4) 

(E5) 


+ J d 3 rd 3 vf (r, v) A (E) (^ j d 3 r'd 3 Vf (r 7 , v') 6 (V 


GM. 1 

- 1 —v ■ 

r' 2 


- n 0 (E) 


which is minimized by requiring 


6S 

Sf 


= - log / — 1 + /3ijj (r; /) + b • L + 2A (E) n (E\ f) = 0. 


(E 6 ) 


(E7) 


12 We correct here two issues in the comparison to the Giirkan & Hopman 
(2007) results made by |Eilon et al.|(2009| Section 4. 3). Fi rst, since j3s and p v 


be estimated by (/3 S ) = y / (/3g(e’ ) ) and ((3 V ) = a/ ((3%). Second, = 


are the rms values, the average values of |Glirkan & Hopman| ^2007) should 


+ rj_, defined by 

Eilon et al. 

2009 1 should be compared to the sum 

+ rj_ of the quantities defined ir 

Giirkan & Hopman 

2007 
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Therefore, the DF that maximizes the entropy is 

/(r,v) = A{E)e 0 ^ +hL , (E8) 

where A(E) = exp[2A(-E)n(-E)], and j3 and b are constants determined by the constraints on n(E ), E+ and L tot . Note that 
an isotropic system must have b = 0 to ensure that the DF does not depend on L, and must also have <3 = 0 since the star-star 
potential depends on L even in an isotropic systenf^] 


E.l. Fluctuation dissipation relation for a spherical symmetric system 

Since stars are assumed to move on Keplerian orbits, it is convenient to work in action-angle coordinates. Choosing the 
^-coordinate in the direction of the total angular momentum, the steady-state is 

p/3-! ti+bL z 

n(E,J,J,)=n(E) ' < E9 > 

For a spherical symmetric system with L to t = 0, the steady-state DF is given by an implicit integral equation 

2 JeW B ’ J ’ n ) 


n (E, J) = n (E) 


J e /3V>(S,J,n)^J2 ■ 


(E10) 


MEP considerations do not require the additional assumptions that go into the FP equation (e.g., diffusion described by a 
Markovian process), and the MEP solution is independent of the path that the system took to reach it from its initial conditions. 
Therefore, the MEP solution must also be satisfied by the FP equation in steady-state, which enforces a connection between the 
DCs, known as the fluct uation-dissipation (F-D) relation. For the symmetries and conserved quantities of the system studied here 
(described by Eq. ( |E10| )), the functional form of the F-D relation is 

2JD je W = E [JDjjeW] . (Ell) 

From this point on we will restrict ourselves to isotropic systems which reach a MEP solution with j3 = 0. This is a solution 
of the form n (E, J) = 2 n (E) J/ 'J%. This means that DCs derived under the assumption of an isotropic background (and fixed 
Keplerian energy) must satisfy the F-D relation 2 JD j = d ( JDjj ) /dJ. In particular, the RR DCs that were used in this study 
obey this relation, as required (|Bar-Or & Alexander|2014|l. 


E.2. Fluctuation-dissipation relation for J-only two-body relaxation 

For NR, the coherence time is shorter than the orbital period (Section|2|, and therefore orbital energy is not conserved. However, 
we argue in Section [3] that the flow pattern in phase-space justifies the approximate treatment of the fast J-diffusion as separate 
from the slower /^-diffusion. We now use this property to show that in that case, the NR J-only DCs also satisfy the fluctuation 
dissipation re latio n, which is a partial test of the validity of the general NR DCs and is used in Section[3] 

Using Eq. ( |C9| > and setting A E = 0, we obtain 

2«Aw|| + Av 2 = 0 . 

Therefore, to first order in Av 2 /v 2 , the local diffusion coefficients are 

The orbit-averaged diffusion coefficients are 

f Va dv 

D xy = {{Ax A y)) 0 =2 {Ax Ay) — , 

J r p v r 

where v r is the radial velocity. Assuming a spherical potential <f> (r), the energy E and angular momentum J are 

E = \ 2 r + J {J 2 /r 2 ) + $ (r), (E16) 

and 

J = rv t , (E17) 

where v t is the transverse velocity. Therefore, 

v r = , (E18) 


(E12) 

(E13) 

(E14) 

(E15) 


13 The L dependence of ^ arises from the eccentricity dependence of the 
enclosed mass seen by the star. This is manifested dynamically by the retro¬ 


grade evolution of the argument of periapse—mass precession. 
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Bar-Or & Alexander 


and 


It then follows that 


and in particular 


d(l/v r ) 11 J 


1 


1 J 


dJ v 2 v r r 2 (r 2 — J 2 /v 2 ) v r v 2 ' 


1 d 


2 JdJ JDjJ 2J - 9 ' (( Aw -l) 




= Dj, 


which therefore proves that the J-only NR DCs indeed satisfy the F-D relation. 


(El 9) 

(E20) 


(E21) 
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